AD-A279  234  ON  PAGE 


form  Approvta 
0MB  No.  0/04  0188 


tr  •  »•>«»  Off  ffiMmo.  incluOMf  (it*  nm*  lof  rr»>*t«M«(  »tiriKtHiit(.  M^rciun*  •..tuivi  a4u  wt«<»t 

vOtlWIHMI  Ol  Ht(>»«lt4(10n  (OflUlWftlt  (Itn  btHOMt  Of  «t»  Olltff  4to*tl  01  ll'H 

•VtMHiion  •<*MaiMft*ft  S*nK*t.  Oi(f<lo«*lt  (Of  •nfm"i«io«  no*f*uon^  onO  Otcwfit  li  iH  itfil«fvin 
4i|«"*»(  luOiiM.  •4o*niraf«  ll«Ouctiaii*fOi*<tig/IM-aiH>.  Wotlwn^lon.  OC  JOiOJ 


1.  AC.  .  ....  . .  3.  REPOBT  TtPt  A 

I  6  May  94  Technical 


4.  TITLE  AND  SUBTITLE 

The  Formulation  of  Quantum  Statistical  Mechanics  Based  on 
the  Feynman  Path  Centroid  Density  IV.  Algorithms  for 
Centroid  Molecular  Dynamles 


6.  AUTH0R(S} 


3.  REPORT  TYPE  ANO  OATES  COVERED 

Technical  5/93-5/94 


S.  FUNDING  NUMBERS 


Jlanshu  Cao  and  Gregory  A. 


7.  PERFORMING  ORGANIZATION  NAME(S)  ANO  AOORESS(E 


Department  of  Chemistry 
University  of  Pennsylvania 
Philadelphia,  PA  19104-6323 


DTIC 


MAY  17 1994 


9.  SPONSORING /MONITORING  AGENCY  NAME(SI  ANO  AOOR£SS{ES) 

( 

Office  of  Naval  Research 
Chemistry  Division 
800  North  Quincy  Street 
Arlington,  VA  22217-5000 


11.  SUPPLEMENTARY  NOTES 


I 

■ 


ONR  N00014-92-J-1243 


PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ONR  Technical  Report 
#9 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


94-14659 


124.  OISTAIBUTION/AVAILABILITY  STATEMENT 

12b.  DISTRIBUTION  CODE 

Approved  for  public  release:  distribution  unlimited 

13.  ABSTRACT  (MdMtmum  200  wor^} 

Numerical  algorittuns  are  developed  for  the  centroid  molecular  dynamics  (cenuoid  MD)  method  to 
calculate  dynamical  time  correlation  functions  for  general  many-body  quantum  systems. 
A{^)roaches  based  on  the  normal  mode  path  integral  molecular  dynamics  and  staging  path  integral 
Monte  Carlo  methods  are  described  to  carry  out  a  direct  calculation  of  the  force  on  the  centroid 
variables  in  the  centroid  MD  algorithm.  A  more  efficient,  but  approximate,  scheme  to  compute  the 
centroid  force  is  devised  which  is  based  on  the  locally  optimized  harmonic  approximation  for  the 
centroid  potentiaL  The  centroid  MD  equations  in  the  latter  mediod  can  be  solv^  with  the  help  of  an 
iterative  procedure  or  through  extended  Lagrangian  dynamics.  A  third  algorithm  introduces  an 
effective  centroid  pseudqx)tential  to  approTdmate  the  fuU  noany-body  Centro^  mean  force  potential 
by  effective  pairwise  centroid  interactions.  Numerical  simulations  for  both  protot]^  models  and 
more  realistic  many-body  systems  are  performed  to  explore  the  feasibility  and  limitations  of  each 
algorithna. 


14.  SUBJECT  TERMS 

t  ' 

Chemical  dynamics;  computer  simulation;  electrochemistry 


IS.  NUMBER  OF  PACES 

53 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


Unclassified 


•(‘..*1  /SAO  O'  .^ao  ssoo 


Unclassified 


Unclassified 


StJnd.rd  form  208  (Rev  2  H’)) 

Ov  ..M'.l  /!'»•  ^ 

iH  I'li 


OFHCE  OF  NAVAL  RESEARCH 

GRANT  N00014-92-J-1243 

R&T  Code  4131065 
Scientific  Offiar:  P.P.  Schmidt 
Technical  Report  No.  9 

The  Formulation  of  Quantum  Statistical  Mechanics  Based  on  the  Feynman  Path 
Centroid  Density.  IV.  Algorithms  for  Centroid  Molecular  Dynamics 

ty 

Jianshu  Cao  and  Gregory  A.  Voth 
Submitted 

to 

Journal  of  Chemical  Physics 

University  of  Pennsylvania 
Department  of  Chemistry 
PhUadelphia,  PA  19104-6323 

May  1994 

Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of  the  United  States 

Government 

This  document  has  been,  approved  for  public  release  and  sale;  its  distribution  is 

unlimited 


(1)* 


Office  of  Naval  Research 
Chemistry  Oivision,  Code  313 
800  North  Quincy  Street 
Arlington,  Virginia  22217-5000 


Defense  Technical  Information 
Center  (2) 

Building  5,  Cameron  Station 
Alexandria,  VA  22314 


Or.  James  S.  Murday  (1) 

Chemistry  Division,  Code  6100 
Naval  Research  Laboratory 
Washington,  D.C.  20375-5000 


Dr.  Robert  Green,  Director  (1) 
Chemistry  Oivision,  Code  385 
Naval  Air  Weapons  Center 
Weapons  Oivision 
China  Uke,  CA  93555-6001 


Dr.  Elek  Lindner  (1) 

Naval  Command,  Control  and 
Ocean  Surveillance  Center 
ROT&E  Division 
San  Diego,  CA  92152-5000 


*  Number  of  copies  to  forward 


Dr.  Richard  W.  Drisko  (1) 
Naval  Civil  Engineering 
Laboratory 
Code  L52 

Port  Hueneme,  CA  93043 


Dr.  Harold  H.  Singerman  (1) 
Naval  Surface  Warfare  Center 
Carderock  Division  Detachment 
Annapolis,  MD  21402-1198 


Dr.  Eugene  C.  Rscher  (1) 
Code  2840 

Naval  Surface  Warfare  Center 
Carderock  Division  Detachment 
Annapolis,  MD  21402-1198 


Dr.  Bernard  E.  Douda  (1) 
Crane  Division 

Naval  Surface  Warfare  Center 
Crane,  Indiana  47522-5000 


Accesion  for 

T 

NTIS  CRA&I 

OTIC  TAB 

□ 

Uiiaiwojnced 

□ 

Justification 

. . . 

By  _ 

Dist:  ibution  f 

Availability  Codes 

Avaif  and  I  or 


Dist  Special 


Ply5 , 

The  Formulation  of  Quantum  Statistical  Mechanics  Based  on 
the  Feynman  Path  Centroid  Density.  IV.  Algorithms  for 
Centroid  Molecular  Dynamics 

Jianshu  Cao  and  Gregory  A.  Voth 
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Abstract 

Numerical  algorithms  are  developed  for  the  centroid  molecular  dynamics 
(centroid  MD)  method  to  calculate  dynamical  time  correlation  functions  for 
generd  many-body  quantum  systems.  Approaches  based  on  the  normal  mode 
path  integral  molecular  dynamics  and  staging  path  integral  Monte  Carlo 
methods  are  described  to  carry  out  a  direct  calculation  of  the  force  on  the 
centroid  variables  in  the  centroid  MD  algorithm.  A  more  efficient,  but  ap¬ 
proximate.  scheme  to  compute  the  centroid  force  is  devised  which  is  based  on 
the  locally  optimized  harmonic  approximation  for  the  centroid  potential.  The 
centroid  MD  equations  in  the  latter  method  can  be  solved  with  the  help  of  an 
iterative  procedure  or  through  extended  Lagrangian  dynamics.  A  third  algo¬ 
rithm  introduces  an  effective  centroid  pseudopotential  to  approximate  the  full 
many-body  centroid  mean  force  potential  by  effective  pairwise  centroid  inter¬ 
actions.  Numerical  simulations  for  both  prototype  models  and  more  realistic 
manv-body  systems  are  performed  to  explore  the  feasibility  and  limitations 
of  each  algorithm. 
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I.  INTRODUCTION 


In  a  series  of  previous  papers,  [1-4]  an  intriguing  perspective  on  quantum  statistical 
mechanics  has  been  extensively  developed.  This  perspective  involves  the  path  centroid  vari¬ 
able  [5]  in  Feynman  path  integration.  [6-12]  (Hereafter,  three  of  the  preceding  papers  [2-4] 
will  be  referred  to  as  Papers  I.  II.  and  III,  respectively.)  Perhaps  the  most  interesting  and 
promising  result  of  the  aforementioned  research  is  an  approximate  method  for  computing 
quantum  d3mamical  time  correlation  functions  called  “centroid  molecular  dynamics"  (cen¬ 
troid  MD).  [1.3.4]  The  central  concept  of  centroid  MD  was  first  introduced  in  Ref.  [1]  and 
then  extensively  analyzed  and  extended  in  Paper  II.  Paper  III.  which  precedes  the  present 
paper,  introduces  a  phase  space  centroid  perspective  and  provides  a  definitive  theoretical 
basis  for  centroid  MD.  Though  a  single  one-dimensional  quantum  particle  has  largely  served 
as  the  numerical  example  for  centroid  MD  in  the  previous  papers.  [1,3,4]  the  method  is  by 
no  means  limited  to  such  a  simple  prototype  model.  As  a  matter  of  fact,  one  prominent 
advantage  of  the  method  is  the  feasibility  of  its  application  to  more  comphcated  physical 
systems  such  as  solvated  quantum  particles,  high  frequency  vibrations  of  condensed  phase 
polyatomic  molecules,  clusters  and  fluids  consisting  of  quantum  particles,  and  nonlinear 
(luantum  phonons.  However,  a  central  challenge  is  to  first  develop  practical  numerical  al¬ 
gorithms  for  computing  centroid  MD  in  general  many-body  quantum  systems.  The  present 
paper  is  devoted  to  this  issue. 

Before  describing  the  algorithms  for  centroid  MD,  it  is  first  necessary  to  briefly  review 
the  key  developments  in  the  earlier  papers.  In  Paper  I.  the  quasiclassical  role  for  the  path 
centroid  variable  and  centroid  density  [5.13-15]  in  defining  the  equilibrium  properties  of 
quantum  systems  was  fully  elucidated.  Then,  building  on  the  insights  of  an  earlier  commu¬ 
nication.  [1]  the  path  centroid  perspective  was  extended  in  Paper  II  to  the  realm  of  quantum 
dynamics.  In  that  paper,  a  centroid-based  perspective  was  uncovered  for  the  calculation  of 
quantum  time  correlation  functions.  A  centroid  time  correlation  function  was  introduced  in 
rhose  papers  as  an  approximation  for  the  Kubo  transformed  [16]  position  correlation  func- 
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tion.  The  centroid  correlation  function  is  thus  related  to  the  usual  quantum  time  correlation 
function  through  the  Fourier  transform  relationship 

C{ij)  =  {hdu /' 2)  [coihi h3uj /2}  +  1]  C'iuj)  1.1) 

where  C*  denotes  the  centroid  time  correlation  function.  [1,3,4)  It  was  argued  that  C*ff) 
should  be  computed  from  classical-like  dynamics  for  the  centroid  variable  on  the  mean 
centroid  potential  energy  surface  such  that 

C’(f)  =  (?c(f)9c(0)>p,  .  (1.2) 

1 

where  qdt)  is  the  centroid  trajectory  for  a  given  particle  degree-of-freedom.  The 
3iV— dimensional  centroid  trajectories  for  a  general  system  are  generated  by  the  effective 
classical  equations  of  motion  [1,3,4] 

,  (1.3) 

where  qdt)  is  the  3iV— dimensional  column  vector  of  centroid  positions,  m  is  the  diagonal 
particle  mass  matrix,  and  Fd^c)  is  the  quantum  mechanical  centroid  mean  force  vector. 
The  latter  quantity  is  expressed  as  the  operation  of  a  3Af— dimensional  gradient  vector  in 
centroid  conhguration  space.  Vc,  on  the  mean  centroid  potential  Vdqc),  i.e.. 

Fdqc)  =  -VcV'c(gc) 

/  •  • .  f  V^r)  6{qc  -  go)  {  W[glO))}  exp{-S{^T)]/h} 

/  •  •  ■  /  P^r)  6(gc  -  go)  exp{-5[^r)]/n} 

Here.  5[^r)]  is  the  imaginary  time  action  functional  [6-9]  and  the  imaginary  time  position 
centroid  vector  defined  in  Eq.  (1.4)  is 

The  effective  temperature-dependent  centroid  potential  Vdqc)  in  Eq.  (1.4)  is  defined  by 

^'c(9c)  =  - A,-sT/n  |pc(gc)/[n{nit/27r/i^3)]^/^|  1.6) 
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where  the  muhidimensionai  position  centroid  density  PcKQc)  is  given  by  [2.5] 

MQc)  =  /•••/  '^^T)^i^c-<lo)exp{~S[q[T)/h]}  (1.7) 

In  the  centroid  correlation  function  [Eq.  (1.2)|,  the  initial  condition  averaging  denoted  by 
(•  •  •)p^  is  formaily  performed  with  the  phase  space  version  of  Eq.  (1.7)  [1,3,4]  in  which  the 
centroid  momentum  distribution  simply  factorizes  into  the  classical  Boltzmann  momentum 
distribution. 

In  Paper  III.  the  centroid  formulation  in  multidimensional  phase  space  allows  one  to 
prove  that  the  centroid  MD  position  correlation  function  is  a  well-defined  approximation 
to  the  Kubo-transformed  position  correlation  function.  [4]  The  agreement  between  the  two 
correlation  functions  can  be  demonstrated  uniformly  to  all  orders  in  a  time  expansion,  with 
an  error  being  proportional  to  the  effective  thermal  width  of  the  particle  and  the  average 
anharmonicity  of  the  mean  centroid  potential.  The  reader  is  referred  to  the  preceding  pa¬ 
per  (Paper  III)  for  an  extensive  discussion  on  this  subject.  In  addition  to  this  analytical 
justification,  an  earlier  justification  [1,3]  was  based  on  the  analytic  continuation  of  correla¬ 
tion  functions  obtained  from  the  effective  harmonic  representation  of  the  centroid  density. 
[2.13-15]  Moreover,  centroid  MD  has  been  tested  numerically  for  some  non- trivial  nonlinear 
systems  and  found  to  give  excellent  agreement  with  the  exact  result  for  the  quantum  posi¬ 
tion  correlation  function.  [1,3]  Several  centroid  MD  strategies  were  also  developed  in  Papers 
II  and  III  to  compute  general  quantxim  time  correlation  functions  of  the  form  (/l(t)B(0)), 
where  A  and  B  are  general  quantum  operators. 

.\lthough  centroid  MD  appears  to  be  a  substantial  breakthrough  in  the  computation 
of  (approximate)  quantum  time  correlation  functions,  the  determination  of  the  centroid 
force  in  Eq.  il.3).  as  defined  by  Eq.  (1.4).  represents  an  algorithmic  challenge  for  realistic 
many- body  s>'srems.  Equation  (1.4)  shows  that  the  centroid  force  is  given  by  the  average  of 
the  local  potential  gradient  over  quantum  path  fluctuations  about  the  constrained  centroid 
variables.  Of  course,  all  numerical  path  integral  simulation  techniques  [17.18]  can  be  adapted 
to  compute  this  average,  but  the  real  issue  is  one  of  computational  efficiency  so  that  the 
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centroid  force  can  be  readily  computed  "on  the  tly"  during  the  real  time  integration  of  the 
centroid  MD  equations  [Eq.  (1.3)|.  This  computation  is  not  trivial  for  quantum  many-body 
systems.  (It  is.  however,  clearly  easier  than  a  ntunericai  "frontal  assault'’  on  the  many-body 
time  dependent  Schrddinger  equation  or  on  real  time  Feynman  path  integrals.  [6]) 

In  the  present  paper,  several  algorithms  will  be  developed  and  explored  for  the  efficient 
computation  of  the  centroid  force  in  the  centroid  MD  equations.  Two  of  the  approaches  will 
be  direct,  or  "brute  force” ,  numerical  path  integral  approaches  [17.18]  in  which  the  path  aver¬ 
aging  in  Eq  1)  is  performed  explicitly.  These  approaches  will  be  based  on  a  normal  mode 
path  integral  MD  algorithm  [19]  and'a  staging  path  integral  Monte  Carlo  (MC)  algorithm. 
[20.21]  with  several  numerical  tricks  introduced  to  speed  up  the  computations.  A  rather 
different,  and  more  efficient,  direction  will  also  be  taken  which  approximates  the  centroid 
force  within  the  ffamework  of  the  locally  optimized  effective  harmonic  perspective.  [2,13-15] 
The  numerical  details  of  solving  the  transcendental  equations  in  the  variational  quadratic 
approximation  are  explored  in  some  detail,  and  two  independent  algorithms  are  developed 
to  optimize  the  variational  parameters  while  simultaneously  propagating  the  centroid  MD 
variables  in  time.  In  addition,  several  schemes  are  proposed  to  simplify  the  evaluation  of 
the  Gaussian  averages  explicit  in  the  wiational  theory.  [2.13-15]  Yet  a  third  algorithm  for 
centroid  MD  involves  the  calculation  of  the  excess  centroid  free  energy  as  a  function  of  the 
separation  between  a  pair  of  particles,  representing  this  free  energy  by  peiirwise  centroid 
pseudopotentials,  and  then  running  centroid  MD  simulations  for  an  effective  many-body 
system  which  interacts  only  through  such  pairwise  interactions.  The  latter  approximation 
reduces  the  full  centroid  mean  potential  hypersurface  to  the  superposition  of  two-body  cen¬ 
troid  potentials  and  thus  ignores  static  many-body  quantum  correlations.  Nevertheless,  the 
latter  approach  will  provide  a  simple  and  effective  alternative  for  centroid  MD  simulations 
of  nearly  classical  or  weakly  interacting  systems.  All  of  the  aforementioned  algorithms  are 
tested  on  several  systems  to  demonstrate  their  overall  feasibility  and  limitations. 

The  present  paper  is  organized  as  follows:  In  Sec.  II.  two  algorithms  are  described  for  the 
direct  numerical  calculation  of  the  centroid  force  during  a  centroid  MD  simulation.  Then. 
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in  Sec.  Ill  an  alternative  and  more  numerically  efficient  method  is  explored  which  is  based 
on  the  eifective  quadratic  (i.e..  approximate)  theory  for  the  centroid  force.  Several  different 
aspects  of  the  implementation  of  this  algorithm  are  also  discussed  in  some  detail.  Section 
IV  then  contains  a  description  of  the  centroid  pseudopotentiai  method,  while  numerical 
implementations  of  all  of  the  algorithms  are  presented  in  Sec.  IV.  Concluding  remarks  are 
given  in  Sec. 


II.  DIRECT  CALCULATION  OF  THE  CENTROID  FORCE 


In  the  discretized  version  of  the  Feynman  path  integral.  [17.18]  the  centroid  force  in  Eq. 
(1.4)  is  written  as 

Ff  )  =  /  •  •  •  /  nUi  (ki  S(gc  -  %}  V'iqi,  •  •  • .  <?p)  exp  [-Sp(qi, . . . ,  qp)/n] 

^  nf=i  dqi  <!»(9c  -  qo)  exp  [Spiqu ....  qp)/n] 

where,  from  the  cyclic  invariance  of  the  trace, 

V'ta . op)  =  (2-2) 

and  the  centroid  \-ariabie  in  discrete  notation  is  given  by 

9o  =  if;*  .  (2.3) 

•*  !=i 

The  discretized  path  integral  action  functional  is  given  by 

Sp(* . 9P)  =  £  ^  (2.4) 

Here,  the  discretization  parameter  is  given  by  P  and.  for  notational  simplicity,  the  expres¬ 
sions  are  written  for  a  single  degree-of-freedom.  An  extension  to  many  degrees-of-freedom 
is  straightforft'ard.  although  notationally  more  cumbersome. 

In  a  path  integral  Monte  Carlo  i  PIMC)  calculation,  the  centroid  force  can  be  readily 
calculated  from  Eq.  (2.1)  by  using  the  importance  sampling  function  exp[—Sp{qi,  —  qpj/^ 
and  pairwise  MC  moves  to  enforce  the  centroid  constraint.  For  a  path  integral  molecular 
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dynamics  (PIMD)  calculation.  (22)  one  defines  fictitious  momenta  p,  tor  each  of  the  quasipar¬ 
ticles  Qi  and  then  runs  an  MD  simulation  with  Hamilton’s  equations  based  on  the  fictitious 
Hamiltonian 


rr  Pi  \  2  ^  I  V  ■ 


where,  in  the  case  of  centroid  MD,  m'  must  equal  m/P  so  that  the  centroid  variable  has  the 


physical  particle  mass  m.  [1,3,41  PIMD  trajectory  adequately  samples  the  canonical 
ensemble,  then  a  particular  time  average  over  that  trajectory  will  yield  the  centroid  force. 


i.e.. 


FciQc)  =  I  dT{V[qiT)])a 

T—oo  T  Jo  P  i^i  dq 


7=?i(t).W=9c 


where  (•  •  •)c  denotes  a  centroid-constrained  path  integral  average  [cf.  Eqs.  (1.4)  and  (2.1)]. 
In  a  PIMD  calculation,  the  centroid  constraint  must  be  enforced  through,  e.g.,  a  holonomic 
constraint. 


While  a  direct  numerical  path  integral  computation  of  the  centroid  force  as  outlined  above 
will  undoubtedly  provide  an  accurate  value  of  the  centroid  potential  surface,  the  relevant 
question  here  is  how  to  affect  such  a  computation  within  the  context  of  the  time-integration 
of  the  centroid  MD  equations  [Eq.  (1.3)].  For  low  dimensional  systems,  the  centroid  force 
can  indeed  be  first  calculated  for  each  point  in  space  and  stored  on  a  grid  in  computer 
memory  to  be  later  recalled  in  the  centroid  MD  calculation.  As  the  dimensionality  of  the 


system  increases,  however,  this  straightforward  procedure  is  no  longer  feasible  due  to  the 
exponential  growth  of  the  required  memory  space.  Consequently,  for  many-body  systems 


one  is  forced  to  repeat  a  PIMD  or  PIMC  calculation  of  the  centroid  force  at  each  centroid 


MD  time  step  (or  over  the  interval  of  each  time  step).  .\s  it  turns  out.  neither  of  the 
straightforward  numerical  path  integral  methods  outhned  above  is  feasible  for  this  purpose. 
In  the  following  two  subsections,  some  more  specialized  algorithms  for  the  calculation  of  the 
force  in  centroid  MD  are  presented. 


A.  Normal  Mode  Path  Integral  Molecular  Dynamics 


One  possible  approach  to  the  computation  of  the  centroid  force  during  a  centroid  MD 
calculation  is  based  on  the  normal  mode  path  integral  molecular  dynamics  (NMPIMD) 
algorithm.  Though  NMPIMD  has  been  discussed  at  length  in  a  different  context.  [19]  it 
\vill  be  presented  here  within  the  framework  of  centroid  MD.  In  the  NMPIMD  algorithm, 
the  centroid  variable  naturally  separates  from  the  other  Feynman  path  modes.  This  feature 
provides  an  important  simplification  of  the  computation  of  the  centroid  force  [Eq.  (1.4)|  and 
in  distinguishing  the  sampling  inherent  in  that  computation  from  the  real  time  propagation 
of  the  centroid  variable.  Moreover,  the  normal  modes  corresponding  to  the  path  fluctuations 
about  the  centroid  can  be  assigned  a  smaller  fictitious  mass  than  the  physical  particle  mass, 
leading  to  a  more  rapid  convergence  of  the  centroid  force. 

The  basic  idea  of  NMPIMD  is  to  diagonalize  some  significant  part  of  the  Feynman  action 
functional  via  a  normal  mode  transformation,  defined  by 

p 

Qj  =  ar,expi-2iTijn/ P)  (2.7) 

n=l 

where  an  is  the  normal  mode  coordinate.  In  order  to  shorten  the  computational  time  to 
perform  the  transformation  in  Eq.  l2.7).  one  can  taJce  P  to  equal  2'^.  where  k  is  an  integer, 
and  then  make  use  of  the  fast  Fourier  transform  technique.  By  virtue  of  the  orthonormality 
relation  cxp[27rim(/  —  n)/P]  =  PSni,  the  above  transformation  will  diagonalize  any 
"quadratic  action  functional.  It  should  be  noted  that  the  normal  mode  propagator  is  different 
from  the  Fourier  path  integral  (FPI)  propagator  [12,23]  which  is  continuous  in  coordinate 
space  but  truncated  in  Fourier  space.  Unlike  the  convergence  of  the  normal  mode  propagator 
which  depends  on  the  discretization  parameter  P.  the  convergence  of  the  FPI  propagator 
depends  on  the  truncation  parameter  of  the  Fourier  modes  kmax- 

Clearly,  the  computational  effort  implicit  in  NMPIMD  relies  on  the  efficiency  of  the  path 
integral  algorithm  which,  in  turn,  depends  on  the  accuracy  of  the  propagator  used  in  the 
sampling.  For  example,  it  is  advantageous  to  include  the  quadratic  part  of  the  potential 
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energy  function  in  the  definition  of  the  NMPIMD  propagator.  The  imaginary  time  action 
functional  for  a  linear  harmonic  oscillator  (LHO)  is  given  in  the  discretized  notation  bv 


oiqi,---qp)/  ^^2h^3sinh(R) 


[{qf  +  qi-i]cosh( R]  -  2g,9,_ii 


where  R  =  hdu.  /  P.  with  u  being  the  LHO  frequency.  If  the  normal  mode  transformation  in 
Eq.  (2. 71  is  appued  to  this  action,  one  obtains 

So(ai . ap)/ft  - 

n=l  -''n 

_  (a«)^  +  (g')^  ^ 

\2  Oli  ‘•>y2  v-yj 


where  A,,  is  the  thermal  width  of  the  n-th  normal  mode,  defined  as 


\  2 

”  2P2[l  -  cos(2;rn/P)]  +  b'i 


(2.10) 


Here,  the  redefined  parameters  according  to  the  exact  LHO  propagator  in  Eq.  (2.8)  are  given 


h^0  sinh{R) 


. ^ 

^  m  R 


(2.11) 


h.  =  hduj 


sinh{R/2) 


(2.12) 


It  is  obvious  from  Eq.  (2.9)  that  the  Gaussian  widths  for  the  real  part  and  imaginary 
part  of  the  first  P/2  —  1  modes  are  given  by  A„.  On  the  other  hand.  ap/i  and  ap  are 
real  variables  with  widths  V2Ap/2  and  \/2Ap,  respectively.  The  normal  mode  variables  from 
(P/2)  +  1  to  P  -  1  can  be  obtained  by  taking  the  complex  conjugate  of  the  first  (P/2)  -  1 
variables  according  to  the  relation  a’  =  ap^n- 

FVom  the  above  equations,  it  becomes  apparent  that  the  normal  mode  ap  in  Eq.  (2.9) 
is  the  centroid  variable  Qq  from  Eq.  (2.3)  which  is  to  be  held  constant  during  a  NMPIMD 
(’alculation  of  the  centroid  force.  One  can  therefore  define  an  effective  ‘'Hamiltonian"  for 
the  remaining  path  modes  by  assigning  them  a  set  of  fictitious  momenta  and  masses,  i.e.. 
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(2.13) 


ffpfe)  =  E 


23Xi)  ‘  P 

where  Pn  is  the  conjugate  fictitious  momentum  for  the  normal  mode  On  and  AV  is  the 
residual  potential,  defined  by 


AV{q)  =  V{q)  -  -rruj^q  -  qcf 


(2.14) 


Molecular  dynamics  can  now  by  computed  from  Hamilton's  equations  based  on  this  effective 
Hamiltonian,  and  the  quasiparticle  coordinates  {gj}  obtained  as  a  function  of  time  from  Eq. 
(2.7).  For  the  sake  of  efficiency,  the  fictitious  masses  are  chosen  as  m„  =  so  that 

all  the  normal  modes  oscillate  with  the  same  unit  frequency.  Otherwise,  the  NMPIMD  time 
step  will  be  determined  by  the  timescale  of  the  fast  mode  and  the  numerical  integration  of 
the  slow  modes  becomes  ineffective.  To  assure  a  canonical  ensemble  average  over  the  normal 
modes,  one  can  attach  a  Nose-Hoover  chain  to  each  of  the  mode  variables.  [24] 

To  our  knowledge,  this  is  the  first  attempt  to  implement  a  NMPIMD  algorithm  in  any 
context.  It  should  be  noted  that  a  very  similar  algorithm,  called  staging  PIMD,  [25]  has 
already  been  implemented  with  considerable  success  for  equilibrium  path  integral  calcula¬ 
tions.  However,  the  NMPIMD  algorithm  is  preferable  for  centroid  MD  calculations  because 
the  centroid  \-ariable  is  naturally  identified  from  the  normal  mode  transformation.  Funher- 
more.  the  programming  of  NMPIMD  is  straightforward.  -A.  comment  is  in  order,  however, 
regarding  the  actual  implementation  of  the  NMPIMD  centroid  force  algorithm  for  multidi¬ 
mensional  systems.  In  such  cases,  the  residual  potential  in  Eq.  (2.14)  should  be  defined  as 
the  deviation  of  the  exact  potential  from  the  quadratic  potential  which  arises  from  only  the 
diagonal  elements  of  the  Hessian  matrix.  The  NMPIMD  expressions  given  herein  will  thus 
basically  be  the  same  for  each  degree  of  freedom.  Otherwise,  a  diagonalization  of  the  Hessian 
raatirx  is  first  required  at  each  centroid  position  in  order  to  implement  the  NMPIMD  algo¬ 
rithm  as  it  has  been  presented.  The  numerical  overhead  inherent  in  such  a  diagonalization  at 
every  centroid  MD  time  step  is  likely  to  be  prohibitive  within  the  context  of  the  dynamical 
calculation.  On  the  other  hand,  although  the  NMPIMD  calculation  of  the  centroid  force 
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must  eventually  converge  at  each  centroid  position,  the  neglect  of  the  off-diagonal  terms  of 
rhe  Hessian  matrix  in  the  definition  of  the  residual  potential  may  slow  the  convergence.  In 
such  cases.  NMPIMD  can  no  longer  be  considered  the  method  of  choice  for  the  computation 
of  the  centroid  force  and  an  alternative  should  be  implemented  (see  the  following  sections!. 

The  NMPIMD  algorithm  can  be  incorporated  into  a  centroid  MD  calculation  in  two  ways; 
The  first  is  to  calculate  the  centroid  force  using  Eq.  (2.6)  at  each  time  step  in  the  centroid 
MD  calculation.  In  this  approach,  the  centroid  force  computation  is  algorithmically  distinct 
from  the  centroid  MD  time  integration.  The  second  approach  is  to.  in  effect,  compute  the 
centroid  force  ‘on  the  fly”  within  the  centroid  MD  algorithm.  [26]  This  calculation  can.  in 
principle,  be  accomplished  by  making  the  fictitious  masses  run  =  w!  jQX\  of  the  normal 
mode  \'ariables  small  enough  so  that  the  centroid  force  is  convergently  computed  on  the 
(longer)  timescale  of  the  centroid  motion.  In  this  case,  the  centroid  force  on  that  timescale 
is  essentially  given  by 


,  1  l^dVig) 


(2.15) 


Pt[  dq 

where  Ate  is  the  time  step  needed  to  accurately  integrate  the  slow  centroid  motion.  Clearly, 
the  average  must  converge  on  the  timescale  Ate-  This  algorithm  is  a  kind  of  “extended 
Lagrangian”  technique  like  the  Car-Parrinello  algorithm.  (27.28j  but  there  is  an  important 
difference.  In  the  latter  case,  the  goal  is  to  have  the  parameters  (e.g.,  plane  wave  coefficients) 
oscillate  quickly  and  “tightly”  arouna  the  minimum  (i.e..  the  adiabatic  ground  state  Bom- 
Oppenheimer  surface).  In  a  centroid  MD  calculation,  however,  the  idea  is  to  have  the  path 
fluctuation  modes  sample  their  full  canonical  equilibrium  distribution  (and  to  do  so  quickly). 
.\gain.  rhe  canonical  sampling  can  be  facilitated  by  attaching  a  Nose -Hoover  chain  to  each 
of  the  path  fluctuation  modes.  [24] 

Several  other  techniques  can  also  be  implemented  to  improve  the  convergence  of  the 
NMPIMD  calculation  of  the  centroid  force  in  a  centroid  MD  calculation.  For  examp ’e. 

(a)  The  reference  potential  for  defining  the  residual  potential  [Eq.  (2.14)]  can  be  updated 
at  each  time  bv  setting  the  effective  harmonic  frequency  [2.13-15]  for  the  centroid  force 
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calculation  to  be  the  frequency  computed  at  the  previous  centroid  time  step.  i.e.. 

.  (2.16) 

The  curv’ature  of  centroid  potential  can  be  evaluated  along  with  the  centroid  force  at  any 
given  centroid  time  step. 

(^b)  If  the  reference  potential  is  constantly  optimized,  the  residual  potential  AV  will  be 
small  and  hence  the  MD  motion  of  the  normal  modes  will  be  dominated  by  the  harmonic  part. 
Therefore^  a  multiple  time  step  procedure  [29]  seems  most  promising  for  such  a  calculation. 
In  such  a  procedure,  the  linear  motion  of  normal  modes  is  integrated  over  a  small  time  step, 
whereas  the  residual  force  is  taken  into  account  at  much  larger  time  intervals. 

(c)  In  the  case  of  an  optimized  reference  potential,  the  effective  anharmonicity  of  the 
residual  potential  should  be  small  and  the  centroid  force  can  be  expanded  linearly  as 

FcK)  =  fete)  +  -  <)c)  (2.17) 

d(ic 

where  the  centroid  force  gradient  dF^Qc) I dqc  is  the  curvature  of  the  centroid  mean  potential, 
given  by 

dFAgc)  ^  d^Kig) 
dqc  dq^ 


The  imagiuar>'  time  integrals  in  this  expression  can  be  readily  expressed  in  discrete  path 
integral  notation  [cf.  Eq.  (2.6)].  All  the  terms  in  Eq.  (2.18)  can  then  be  evaluated  along 
with  the  centroid  force  for  a  given  centroid  position.  Prom  Eq.  (2.18),  the  deviation  of  the 
centroid  force  F,:{g'c)  at  a  new  configuration  q'  from  the  exact  centroid  force  Fdqc)  at  the  old 
configuration  Oc  can  be  assigned  some  tolerance,  above  which  a  new  NMPIMD  is  initiated 
to  calculate  a  new  centroid  force,  etc.,  for  the  new  configuration.  In  addition,  the  predicted 
(‘entroid  force  FAq't.)  could  also  be  compared  to  the  classical  force  at  <7',  the  latter  being  used 
if  the  difference  from  the  predicted  centroid  force  falls  below  some  tolerance. 

The  above  list  of  numerical  “tricks’’  is  by  no  means  exhaustive.  There  is  a  great  deal  of 
room  for  future  algorithmic  development  in  centroid  MD. 
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B.  Staging  Path  Integral  Monte  Carlo 


As  was  pointed  out  earlier,  the  centroid  force  can  be  evaluated  by  either  a  PIMC  or  PIMD 
algorithm.  A  major  drawback  of  the  NMPIMD  method  is  the  consumptive  computational 
effort  involved  in  the  normal  mode  transformation  (Eq.  (2.7)!  at  each  time  step  to  transform 
back  to  the  quasiparticie  coordinates.  This  transformation  can  become  rather  problematical 


if  the  number  of  path  integral  discretizations  P  is  large.  As  an  alternative,  a  centroid  MD 
algorithm  based  on  staging  PIMC  (20,21)  is  explored  in  the  present  subsection  which  is 
particularly  advantageous  if  the  number  of  discretizations  is  large. 

The  particular  implementation  of  the  staging  algorithm  in  the  present  work  follows  closely 
the  method  discussed  at  length  in  Ref.  [30j.  However,  it  is  helpful  to  first  incorporate  the 
effective  harmonic  reference  system  into  the  definition  of  the  transition  probability  distri¬ 
bution  ftmction  for  selecting  trial  configmations  of  the  path  integral  quasiparticie  chain. 
(20,21,30)  A  LHO  thermal  propagator  is  given  by 


^0(9, 9',  r)  = 


(  m  R 


m  R 
2h^0  sinh{R) 


[{q'^  +  fyosh{R) 


with  R  =  riwT.  It  can  be  shown  that  the  transition  probability  function  for  a  point  q 
intermediate  between  two  points  qi  and  92  is  given  by 


Piq) 


Pi){qi,q,Ti)po(q,q2,  T2) 
Po(9i,92,n  H-'rz) 


1 


2A2 


(2.20) 


where  q'  is  the  center  of  the  Gaussian  distribution,  given  by 


,  _  qisinh(R2)  +  q2sinh(Ri) 

^  nnhiRi  -f-  R2) 

and  A  is  the  width  of  the  Gaussian  distribution,  given  by 

1  _  mu  [1  1 

A2  h  i  tanhiRi)'^  tanh{R2)\ 


(2.21) 


(2.22) 
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The  above  expressions  can  be  reduced  to  the  more  standard  free  particle  case  [30]  in  the 
limit  oi  R  —  0.  In  the  context  of  the  staging  algorithm.  Eq.  (2.20)  defines  Levi's  recursive 
scheme  for  directly  sampling  conditional  Brownian  motion  paths  which  are  then  subjected 
to  the  acceptance  probability,  determined  by  the  residual  probability 


P(gi>g2.-  •  •  -Qp) 
Po(9i.92.----9p) 


(2.23) 


where  AV  is  given  by  Eq.  (2.14). 

An  algorithmic  challenge  is  to  “smoothly”  incorporate  staging  PIMD  into  the  centroid 
MD  time  integration.  An  accurate  centroid  MD  calculation  may  require  a  small  enough 
time  step  so  that  the  force  difference  in  the  numerical  MD  integrator  becomes  similar  to 
the  statistical  error  in  the  PIMC  sampling  of  the  centroid  force.  To  deal  with  this  issue,  a 
single  centroid  MD  time  step  At  can  be  divided  into  N^c  segments.  As  a  new  MC  chain 
configuration  is  sampled  for  a  fixed  centroid,  the  centroid  can  then  be  moved  according 
to  a  much  smaller  time  step  6t  —  using  the  force  of  the  instantaneous  MC  path 

integral  configuration.  The  Nmc  staging  PIMC  samplings  which  are,  in  principle,  to  be 
performed  at  the  beginning  of  each  time  step  At.  now  become  evenly  distributed  on  a  fine 
MD  grid  of  spacing  6t  so  that  the  PIMD  evolves  *‘on  the  fly”  with  the  centroid  MD  motion. 
Such  an  approach  should  be  more  numerically  effective  than  simply  computing  the  centroid 
force  at  the  beginning  of  each  time  step.  Though  the  fluctuations  of  the  centroid  force 
can  be  obser\'ed  on  the  fine  timescale  6t.  the  centroid  trajectory  becomes  smooth  on  the 
larger  timescale  Af.  To  some  degree,  this  behavior  is  analogous  to  the  stochastic  motion 
of  a  heavy  particle  solvated  in  a  light  solvent,  namely,  Brownian  motion.  [31]  The  manifest 
diffusion  of  a  Brownian  particle  consists  of  numerous  collisions  with  light  solvent  particles. 
.\s  the  Brownian  particle  becomes  heavier,  the  jigsaw  behavior  becomes  less  detectable. 
This  is  precisely  the  behavior  of  the  centroid  variable  in  centroid  .MD.  In  the  extreme  limit 


of  Nmc  —  oc.  one  recovers  the  mean  force  average  according  to  the  PIMD  case  of  T  —  xi 
in  Eq.  (2.6).  In  this  case,  the  force  fluctuations  on  the  centroid  will  become  completely 
smoothed  out  so  that  the  centroid  feels  the  mean  centroid  MD  force. 
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The  complete  staging  PIMC  centroid  MD  algorithm  can  be  summarized  as  follows: 

(1)  Generate  a  MC  chain  segment  (i.e..  moves  of  a  subset  of  quasiparticles)  by  the  staging 
method  based  on  Eq.  (2.20). 

(2)  Uniformly  adjust  the  positions  of  all  quasiparticles  to  fix  the  centroid  position. 

(3)  Accept  or  reject  the  new  configuration  by  comparing  the  change  in  the  residual  potential 
AU  according  to  Eq.  (2.23). 

(4)  Calculate  the  instantaneous  centroid  force  for  the  resulting  quasiparticle  configuration 
and  move  the  centroid  accordingly  for  a  small  time  step  6t  =  i^t/Nssc  in  the  velocity  Verlet 
algorithm. 

(5)  Repeat  steps  ( l)-(4)  Nmc  times  to  complete  one  centroid  MD  time  step  At  within  the 
framework  of  the  centroid  MD  integration. 

The  complete  centroid  MD  simulation  should  to  be  carried  out  for  different  values  of 
jV.v/c  so  that  convergence  of  the  the  mean  centroid  force  can  be  assured.  Moreover,  in  the 
staging  part  of  the  algorithm  [part  (1)|,  the  number  of  quasiparticles  moved  in  each  chain 
segment  should  be  adjusted  to  achieve  a  reasonable  MC  acceptance  ration  in  part  (3).  As 
in  the  NMPIMD  method,  the  effective  harmonic  reference  propagator  in  Eq.  (2.19)  for  a 
multidimensional  system  can  be  defined  from  the  diagonal  terms  of  the  Hessian  matrix  in 
(jrder  to  avoid  a  costly  diagonalization  procedure. 


C.  General  Considerations  for  the  Direct  Calculation  of  the  Centroid  Force 


Clearly,  the  path  integral  average  of  the  centroid  force  at  each  time  step  is  the  major 
workload  of  both  path  integral  algorithms.  Forttmately,  the  effort  to  sample  the  centroid 
force  is  much  reduced  in  regions  of  low  anharmonicity,  for  nearly  classical  systems,  and/or 
for  weakly  interacting  particles.  This  assertion  can  be  demonstrated  by  explicitly  expressing 
the  centroid  force  fluctuations  as 


= 


_  1  .  f«y2 

23m  dui 


(2.24) 
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where  f{q)  =  V''(g).  In  the  above  equation,  the  first  equaiiiy  comes  from  the  rearrangement 
of  Eq.  C--18)  and  the  second  expression  arises  from  the  locally  optimized  effective  harmonic 
approximation.  ,2.13-15]  The  resulting  expression  indicates  that  the  quantum  fluctuations 
of  the  centroid  force  are  proportional  to  (/")c.  which  arises  from  the  averaged  ajxhaxmomcity 
of  the  potential,  and  from  Cc  which  measures  the  effective  quantum  thermal  dispersion  of 
the  particle.  [2-4] 

As  a  final  point  in  this  section,  it  is  noted  that  in  a  direct  path  integral  centroid  MD 
calculation  one  can  also  explicitly  compute  the  centroid-constrained  path  integral  average  of 
a  quantum  operator  at  each  time  t  in  the  time  integration  (i.e..  a  Gaussian  approximation  for 
the  averaged  operator  [3,4.32]  can  be  avoided).  Such  averages  are  required  in  the  method 
called  "centroid  MD  with  semiclassical  operators"  developed  for  the  computation  of  the 
approximate  centroid  correlation  function  [cf.  Eqs.  (4.10)  and  (4.11)  of  Ref.  [4]],  so 

it  is  fortunate  they  can  be  obtained  so  directly. 

III.  EFFECTIVE  HARMONIC  COMPUTATION  OF  THE  CENTROID  FORCE 

.\n  alternative,  highly  eflacient.  algorithm  for  computing  the  centroid  force  in  centroid 
MD  makes  use  of  the  variational  effective  quadratic  theory  for  the  centroid  potential. 
[1-4.13-15]  This  approach,  which  is  approximate,  represents  the  centroid  potential  as  a 
locally  optimized  quadratic  potential  centered  at  the  centroid  position.  The  effective  poten- 
- —  tial  is  then  constantly  updated  as  the  centroid  propagates  in  the  MD  time  evolution. 

Following  the  multidimensional  formalism  presented  in  Paper  III,  an  effective  quadratic 
centroid  force  constant  matrix  can  be  introduced  as 

Kr(qc)  =  'K(qc  +  |))c 

=  . .  ^  f  d^K(qc  +  ^)exp[-i}-C~^(0.qc)-^/2\  (3.1) 

Ydet[27rC^(0,9c)]  '• 

where  the  Gaussian  width  factor  matrix  Cc(0,^),  in  this  case,  is  the  position  sub-block  of 
'he  generalized  centroid-constrained  correlation  function  matrix  in  Eq.  (2.5)  of  Paper  III. 
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(3.2) 


The  Hessian  matrix  in  Eq.  (3.1)  is  defined  by 

K((7;  =  -T-T—  • 

dQidqj 

On  the  level  of  the  effective  quadratic  approximation,  the  Gaussian  width  matrix  is  formally 


expressed  as. 


C,(0,9d  =  5;[;3mnJ  +  ,aK,te)] 


where  m  is  the  3iV— dimensional  particle  mass  matrix  and  Qn  =  27m/ A/?.  A  centroid- 
dependent  unitary  matrix  U(^)  can  be  found  to  diagonalize  the  mass-scaled  centroid  force 
constant  matrix  Ke(^),  giving  the  eigenfrequencies 

Un^:)K.(ft:)U(9,)  =  1.5  .  (3.4) 

where  jjI  is  the  column  vector  of  centroid-dependent  eigenvalues  and  I  is  the 
37V— dimensional  Identity  matrix.  The  Gaussian  width  factor  matrix  in  Eq.  (3.1)  can  be 
determined  from  the  relation 

c.(0,f.)  =  Ute)  [l-S®]  O'te)  .  (3.5) 

where  the  indi\idual  elements  of  the  normal  mode  thermal  width  factor  vector  are  given  by 

^  I  tanh(h0uJc.i/2)  / 

and  the  mass-scaled  unitary  matrix  U{qc)  is  given  by 

.  (3.7) 

Thus,  the  set  of  optimized  frequencies  are  variationally  obtained  as  the  self-consistent 
solution  to  the  transcendental  matrix  equations  (3.1) -(3.8)  for  a  given  centroid  position. 


A.  Evaluation  of  the  Gaussian  Averages 


A  central  algorithmic  challenge  in  using  the  effective  harmonic  theory  is  the  computation 
of  the  Gaussian  averages  inherent  in  the  variational  expression  for  the  centroid  potential. 
[2-4.13-15]  The  latter  expression  is  given  by 
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-•  1  /*  — 

{V{qc  +  ^))c  =  d^Viqc  +  q)exp\-fj-C:^{Q.qc)-q/2]  (3.8) 

ydet(27rC,(0,9c)r  ^  ^ 

From  this  expression,  one  can  obtain  the  etfective  harmonic  approximation  to  the  centroid 
force,  given  by 


Wc)  =  -  V  e(V'(q,  +  I)), 


(3.9) 


where  the  notation  "•  •  •  Icc”  indicates  that  the  centroid  derivatives  in  the  gradient  do  not 
act  on  the  width  matrix  Cc-  The  elements  of  the  centroid-dependent  force  constant  matrix 
in  Eq.  (3.1)  can  be  similarly  expressed  as 

(K.(4;)ly  =  +  (310) 

It  should  be  pointed  out  that  the  Gaussian-averaged  potential  at  the  centroid  qc  in  Eq. 
(3.8)  is  different  from  the  effective  centroid  potential  introduced  in  Eq.  (1.6).  However,  the 
centroid  force  in  Eq.  (3.9)  can  be  shown  to  be  the  effective  harmonic  approximation  to  the 
exact  centroid  force  in  Eq.  (1.4). 

Many  systems  can  be  well  described  by  pairwise  (or  site-site)  potentials.  If  such  is  the 
case,  the  Gaussian  average  in  Eq.  (3.8)  can  be  expressed  in  terms  of  the  summation  over  all 
pair  interactions,  i.e.. 

{V{qc  +  Q))c  =  •  (3-11) 

i=i  j>t 

where  is  the  three-dimensional  vector  connecting  the  z-th  and  j-th  particles.  For  a  specific 
pair  of  interactions,  the  other  degrees  of  freedom  can  be  integrated  out  of  the  average,  leading 
to  a  Gaussian  average  in  a  lower-dimensional  space: 


{v(fij))c  = 


y'det(27rC;,j(gc)] 


^  f  df  f)  exp  [-r  .  •  r/2]  .  (3.12) 

7^1  J  ^ 


where  the  3-dimensional,  centroid-dependent  submatrix  ,j(^;)  can  be  reduced  from  the  full 
matrix  width  factor  Cc,,j(0,^).  The  expression  for  C'y(^c)  and  its  derivation  can  be  found 
in  the  Appendix.  The  above  simplification  reduces  the  computational  effort  considerably  for 
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systems  described  by  p€dr  potentials.  It  must  be  noted,  however,  that  the  submatrix  [Qc) 
still  depends  on  the  global  configiuration  of  the  particle  centroids.  In  practice,  however,  it 
can  be  assumed  that  the  effective  widths  of  a  pair  of  particles  will  not  be  affected  by  a 
change  in  the  position  of  those  particles  if  they  are  beyond  some  cut-off  separation. 

The  Gaussian  average  in  Eq.  (3.8)  can  only  be  carried  out  analytically  for  polynomials. 
Gaussian  potentials,  exponentials,  and  their  combinations.  Gaussian  averages  of  other  func¬ 
tions  cannot  be  expressed  in  closed  analytical  forms.  The  numerical  integration  of  even  a 
three-dimensional  Gaussian  average  [cf.  Eq.  (3.12)  is  time  consuming  and  would  reduce  the 
efficiency  of  the  effective  harmonic  algorithm  for  centroid  MD.  These  difficulties  are  magni¬ 
fied  when  several  iterations  are  required  to  achieve  a  convergent  solution  to  the  variational 
transcendental  equations.  The  best  strategy  therefore,  is  to  represent  the  physical  potential 
by  a  set  of  functions  which  can  be  Gaussian  averaged  analytically.  Several  such  approaches 
are  disctissed  below  which  may  be  appropriate  under  different  circumstances. 

1.  Taylor  Series  Expansiort 

In  the  nearly  classical  regime,  the  Gaussian  width  matrix,  i.e..  the  centroid-constrained 
propagator  matrix  Cc{0,  (^),  is  a  small  quantity  so  the  Taylor  expansion  of  the  force  constant 
matrix  K(^  -f  q)  in  Eq.  (3.1)  can  be  truncated  at  the  first  nonvanishing  term,  giving 

K.®)  Kte)  +  ic.(0,fJ:(9,a,K®l  .  (3.13) 

where  didjK{qc)  is  a  fourth-order  tensor.  The  centroid  force  expression  is  of  a  similar  form, 
but  with  the  matrices  Kc  and  K  replaced  by  the  vectors  Fc  and  F.  respectively.  For  pairwise 
interactions,  the  force  constant  expression  can  be  reduced  to  one  in  three-dimensions  v.ith 
C'  replacing  C.  in  Eq.  (3.13). 
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■2.  Gausstan  Fit  of  the  Potential 


Few  physical  potentials  resemble  a  Gaussian  function.  However,  by  least-squares  ritting 
rhe  potential  data.  Gaussian  functions  can  be  used  to  represent  any  physical  potential,  i.e.. 


^’(9)  =  exp[-(9-d,)- A,  ‘  •  (v-d,)/2J 


V3.14) 


where  A;  is  the  Gaussian  width  matrix  for  a  given  index  1.  Since  the  set  of  Gaussians  does  not 
rigorouslv  form  a  complete  set.  Eq.  (3.14)  is  not  an  expansion  but  is  instead  an  approximate 
fitting  of  representative  points  in  the  range  important  for  the  particle  interactions. 

With  the  above  e.xpression  in  hand,  one  can  easily  complete  the  Gaussian  integration  in 


Eq.  (3.8).  giving  e.xplicitly 


exp  [- (qc  -  <ff)  •  ( A/ -H  Cc)  ^•(^:-t^)/2j  (3.15) 


which  indicates  that  the  Gaussian  averaging  broadens  the  original  Gaussian  fitting  width.  If 
the  interaction  is  described  by  a  pairwise  central  force,  the  Gaussian  average  can  be  carried 


out  for  each  pair  in  Eq.  (3.11),  giving 

(t'(rq))c  =  det^f  ^  1)  ‘  ’ 


(3.16) 


Here,  a?  =  U(I.  where  ai  is  a  constant  for  a  pair  potential.  In  particular,  the  pair  potential 
between  two  sites  is  fit  bv  Gaussians  as 


Ar)  =  £  7/  e 


(3.17) 


This  e.xpression  represents  the  simplest  possible  Gaiissian  fit  in  which  the  origin  of  each 
Gaussian  is  taken  to  be  r  =  0.  The  fitting  procedure  can  obviously  be  generalized  to  use 
distributed  Gaussians. 


3.  Hermtte  Polynomial  Representation  of  Pair  Potentials 


In  principle,  any  pair  potential  can  be  rigorously  expanded  in  terms  of  the  complete 
■•et  of  eigenfunctions  of  the  linear  harmonic  oscillator.  These  functions  are  the  product  of 
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Hermite  polynomials  and  a  Gaussian  function.  Thereby,  the  Gaussian  averages  of  those  basis 
functions  in  Eq.  (3.8)  can  be  written  as  combinations  of  Gaussian  functions,  polynomials, 
and  error  functions.  This  approach,  because  of  its  rigor,  has  an  ad\’antage  over  the  previous 
two  approximations.  However,  it  will  be  more  complicated  both  analytically  and  numerically. 
The  details  of  this  approach  will  be  left  for  future  investigations. 

Plane  Wave  Representation  of  Pair  Potentials 

A  given  pairwise  interaction  between  particles  or  molecular  sites  can  also  be  expanded 
in  terms  of  a  discrete  set  of  planewaves.  In  this  case,  the  pair  terms  in  Eq.  (3.11)  can  be 
expressed  as 

.  (3.18) 

k 

where  is  the  three-dimensional  Fourier  transform  of  the  pair  potential  for  the  discrete 
wavevector  k.  .\fter  performing  the  Gaussian  averaging  explicit  in  Eq.  (3.12),  one  obtains 
for  the  terms  in  Eq.  (3.11)  the  expression 

{v{fij))c  =  ^  V'^eip  [tfc  •  fc.<j  -  k-  C'(qc)  •  k/2\ .  (3.19) 

k 

This  approach  is  both  mathematically  rigorous  and  simple  to  implement,  but  it  can  become 
computationally  expensive  if  a  large  number  of  planewave  basis  functions  are  required  to  fit 
— the  pair  potential. 

Having  discussed  strategies  for  computing  the  Gaussian  averages  explicit  in  the  effec¬ 
tive  harmonic  theory,  two  algorithms  will  now  be  described  for  carrying  out  centroid  MD 
calculations  using  the  effective  harmonic  representation  of  the  centroid  force. 

B.  An  Extended  Lagrangian  Method 

The  optimal  parameters  for  the  effective  harmonic  approximation  to  the  centroid  poten¬ 
tial  Vriqc)  [EQ-  ( 1-6)1  ^*^6  determined  by  a  variational  minimization  of  Vrfqc)  w'ith  respect 
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to  those  parameters.  [2.13-15]  This  mmimization  is  accomplished  with  the  help  of  a  Gibbs- 
Bogoliubov  variational  principle  for  the  centroid  potential  (i.e..  centroid  free  energj^j.  ;33| 
Therefore,  an  extended  Lagrangian  technique  [27.28]  can  be  employed  to  minimize  the  eifec- 
tive  centroid  potential  and  the  resulting  centroid  force  while  simultaneously  propagating  the 
centroid  variables  in  time.  The  similarity  between  centroid  MD  with  the  effective  harmonic 
theory  for  the  centroid  force  and  the  Car-Parrineiio  (CP)  extended  Lagrangian  method 
[27,28]  is  clear:  The  centroid  motion  is  the  analogue  to  the  CP  nuclear  motion,  the  varia¬ 
tional  parameters  in  the  effective  harmonic  theory  are  the  analogues  to  the  CP  planewave 
coefficients,  and  the  centroid  potential  is  similar  to  the  ground  state  CP  adiabatic  energy 
fimctionai.  In  order  to  outline  the  extended  Lagrangian  method  as  it  applies  to  centroid 
MD.  the  following  discussion  will  be  restricted  to  a  one-dimensional  problem  for  simphcity. 
The  algorithm  for  multidimensional  systems  is  described  at  the  end  of  the  section. 

In  the  effective  harmonic  approximation,  [2,13-15]  the  centroid  potential  in  Eq.  (1.6)  is 
given  by 

=  -fcBrin[(6/2)/sinh(6/2)|  +  (V"(9c  +  g))c  -  mu;2Cc(0,gc)/2  .  (3.20) 

where  b  =  hdue,  with  Uc  being  the  frequency  of  the  effective  harmonic  potential,  and 
the  effective  harmonic  expression  for  Cc(O.gc)  is  given  in  Eq.  (3.31.  The  parameter  ^  is  the 
v'ariational  parameter  which  is  optimized  to  minimize  the  centroid  potential.  The  parameter 
can  be  chosen  to  be  the  centroid  frequency  ujc  or  the  width  factor  Cc. 

At  this  point,  one  can  introduce  the  extended  Lagrangian 

+  -  Vciqc,0  ■  (3-21) 

where  ^  is  the  velocity  associated  with  the  variational  parameter  and  is  its  fictitious 
mass.  The  above  Lagrangian  forms  the  basis  for  the  centroid  MD  equations  in  which  the 
variational  parameter  ^  and.  in  turn,  the  centroid  force  are  computed  simultaneously  ^'^th 
the  centroid  motion. 

If  the  variational  parameter  ^  is  taken  to  be  the  frequency  lu'^.  one  obtains  the  nonlinear 
equation  for  the  hctitious  force  on  that  variable,  i.e.. 
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(3.221 


awc 

=  -  [{d'iViq^  -h  q))c  -  mo;;] 
muj^ 

where  dc  is  the  partial  derivative  with  respect  to  the  centroid  variable,  and 

g{b)  =  - - - - \ -  . 

b  4  tanh(6/2)  8  sinh‘^(6/2) 


(3.23) 


The  centroid  constrained  average  in  Eq.  (3.22),  which  depends  on  Uc,  is  given  from  Eq. 
(3.1).  The  variational  solution  for  the  optimized  centroid  free  energy  [2,13-15]  can  obtained 
by  setting  Eq.  (3.22)  equal  to  zero  and  self-consistently  solving  for  the  variational  parameter 
at  a  given  centroid  position.  This  variational  solution  provides  the  initial  condition  for  the 
^  trajectory  generated  by  the  Euler-Lagrange  equations  based  on  the  extended  Lagrangian 
[Eq.  (3.21)]. 

If  the  variable  ^  in  the  extended  Lagrangian  is  instead  taken  to  be  the  Gaussian  width 
factor  Cc,  one  ootains  the  fictitious  force 


^JeKQci^c)  =  "qq 


(3.24) 


in  which  the  following  relation  has  been  used; 


duj'^  mwl 


(3.25) 


Since  the  width  factor  Cr.  is  in  this  case  the  dynamical  variable  instead  of  Uc,  the  centroid 
frequency  can  be  obtained  from  Q  through  the  one-diraensionai  versions  of  Eqs.  (3.5)  and 
(3.6).  As  in  the  case  of  ^  =  cjc,  the  optimal  variational  value  for  Cc  is  obtained  by  setting 
Eq.  (3.24)  equal  to  zero  and  self-consistently  solving  for  Cc-  For  a  given  centroid  initial 
condition,  this  value  of  Cc  would  provide  the  initied  condition  for  the  variable  Cc  in  the 
extended  Lagrangian  simulation. 

Though  the  choice  of  the  variational  parameter  ^  in  Eq.  (3.21)  is  equivalent  for  the 
minimization  of  the  centroid  free  energy,  Eq.  (3.24)  is  much  easier  to  generalize  to  multidi¬ 
mensional  systems,  giving 


23 


(3.26) 


dCc 


where  the  relevant  quantities  here  are  defined  in  Eqs.  (3.1)— (3.8).  For  a  system  of  3A' 
degrees  of  freedom,  there  axe  3A'’(3jV  -}- 1  )/2  independent  variational  variables  (i.e..  elements 
of  the  matrix  Cc).  The  multidimensional  trajectory  from  the  extended  Lagrangian  generates 
the  width  factor  matrix  Cc{t)  as  a  function  of  time.  In  order  to  compute  the  force  in  Eq. 
(3.26),  the  matrix  trajectory  Cdt)  is  used  in  Eq.  (3.1)  to  compute  the  first  term  on  the 
right-hand'Side  of  Eq.  (3.26).  The  second  term  is  extracted  from  Ce{t)  by  inverting  Eq. 
(3.5)  to  get  the  elements  oiiiqc)  and-  in.  turn,  by  solving  for  the  vector  from  Eq.  (3.6). 
The  latter  vector,  along  with  the  transformation  matrix  U((fc),  is  then  used  in  the  second 
term  on  the  right-hand-side  of  Eq.  (3.26). 

The  essence  of  the  extended  Lagrangian  method  [27,28]  is  to  cLoose  such  small  fictitious 
masses  {m^}  in  the  multidimensional  version  of  Eq.  (3.21)  that  the  fictitious  variables  {^} 
rapidly  oscillate  around  the  minimum  of  the  centroid  free  energy  surface  Eq,  (3.20).  The 
centroid  variables  should  then  exhibit  an  adiabatic,  conservative  motion  because  there  is 
little  energy  exchange  between  them  and  the  fictitious  degrees  of  freedom.  One  useful 
technique  is  to  associate  a  Nose -Hoover  chain  [24]  to  each  fictitious  variable  to  keep  it  from 
heating  up. 


C.  An  Iterative  Method  for  the  Centroid  Force 

While  the  e.xtended  Lagrangian  method  discussed  above  is  quite  elegant,  it  turns  out 
that  a  simple  iterative  solution  is  quite  feasible  for  the  multidimensional  effective  harmonic 
variational  equations  [i.e..  Eq.  (3.26)  set  equal  to  zero].  This  algorithm  is  rather  straight¬ 
forward:  Solve  the  transcendental  equation  iteratively  at  a  fixed  centroid  until  convergence, 
calculate  the  centroid  force  and  move  the  centroid  accordingly  in  the  MD  integrator,  and 
repeat  the  procedure  for  each  time  step.  .Many  fewer  iterations  are  required  if  the  optimal 
[parameters  from  the  previous  time  step  are  used  as  the  initial  guess  at  the  current  time  step. 
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Provided  the  centroid  displacement  is  small,  the  convergence  of  the  iterative  scheme  should 
be  fast. 

One  difficulty  associated  with  the  iterative  method  occurs  when  negative  curvatures  are 
present  in  the  classical  Hessian  matrix  for  a  given  centroid  configuration  (i.e..  at  inflection 
points  or  in  a  barrier  region).  If  the  initial  guess  of  parameters  is  poor,  the  iterative  scheme 
may  not  converge,  even  though  the  Gibbs-Bogoliubov  variational  principle  [33]  insures  the 
centroid  potential  is  ffiaite  and  therefore  that  a  solution  must  exist.  If  such  a  problem  arises 
in  the  convergence  of  the  iterative  algorithm  at  a  given  time  step,  the  iterative  process  can 
be  re-started  unth  a  smaller  value  of  h  and  that  parameter  can  then  be  adjusted  until  the 
neighborhood  of  a  fixed  point  is  located.  For  h  =  0.  the  classical  force  is  obtained  which, 
of  course,  is  well-defined.  It  should  be  noted  that  the  effective  quadratic  theory  is  being 
implemented  here  within  the  context  of  the  local  centroid  force  for  each  time  step  in  the 
centroid  MD  integration.  This  local  instantaneous  centroid  potential  is  not  being  used  to 
extrapolate  the  dynamics  to  long  times.  If  such  were  the  case,  unphysical  instabilities  in  the 
centroid  dynamics  might  indeed  be  present. 

IV.  A  PSEUDOPOTENTIAL  METHOD  FOR  THE  CENTROID  FORCE 

For  a  system  consisting  of  particles  interacting  through  pairwise  potentials,  the  three- 
body  quantum  correlations  will  diminish  as  the  distances  increase.  Therefore,  at  relatively 
dilute  densities  only  the  pairwise  interactions  need  to  be  treated  quantum  mechanically. 
This  approximation  in  the  centroid  theory  is  analogous  to  a  treatment  of  many-body  elec¬ 
tronic  polarization  interactions  through  only  a  pairwise  London-like  attraction  force,  thus 
neglecting  the  higher-order  many-body  dispersion  terms.  The  analogy  in  this  case  is  drawn 
between  the  quantum  path  fluctuations  about  the  particle  centroids  and  the  electronic  fluc¬ 
tuations  about  the  mean  particle  wavefunctions.  The  latter  many-body  dispersion  problem 
has  been  systematically  studied  [34]  through  a  quantum  Drude  oscillator  model  [35]  and  it 
was  found  that  at  low  densities  the  two- body  interaction  is  by  far  the  dominant  one. 
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The  exact  excess  free  energy  of  the  centroids  of  a  many  body  system  can  be  approximated 
as  the  superposition  of  the  excess  free  energy  of  the  pair  centroid  interactions,  i.e.. 

.V  .V 

^ciQc)  ^  •  (4.1) 

«=l  j>i 

where  the  effective  centroid  pair  potential  Veff{rc.tj)  is  computed  from  the  expression  for 
two  particles: 

exp[-3veff{r ^,12)] 

/  •  •  •  /  Wi('r)'P92(r)  6(rc.i2  -  Ha)  eip  {-5[gi(r),  q2{T)\/h} 

/•••/  Wi(r)W2('r)6(rc.i2-n2)eip{-5/p[9i(r),g2(T)]/n} 

Here.  q\  and  q2  are  the  three-dimensional  vectors  for  par  :  xes  T”  and  “2”.  respectively. 

«5[9i(’')-92(^)]  is  the  action  functional  for  the  interacting  .;articles.  Sfp[qi{T),q2{T)]  is  the 

free  particle  action  fimctional,  and 

fi2  =  dr  \qi{T)  -  q2{T)\  .  (4.3) 

The  effective  centroid  pair  potential  Ve}f{rc,ij),  which  is  a  central  force  and  a  function  of 
temperature  as  well  as  mass,  can  be  calculated  accurately  by  the  direct  PIMC  sampling 
methods  discussed  in  Sec.  II  and  saved  on  a  one-dimensional  grid.  The  centroid  MD  equa¬ 
tions  [Eq.  (1.3)j  can  then  be  integrated  using  the  centroid  force  computed  from  the  pairwise 
centroid  pseudopotential.  Obviously,  this  scheme  will  be  most  effective  for  nearly  classical 
systems  where  many-body  quantum  correlations  are  negligible,  or  for  weakly  interacting 
systems  where  there  is  a  low  probability  for  three  particles  to  be  mutually  within  the  some 
relevant  interaction  range.  For  highly  condensed  phases  of  strongly  quantized  particles,  the 
compact  structure  insures  that  there  will  be  several  neighbors  for  each  quantum  particle 
and  thus  the  contribution  from  many-body  quantum  correlations  becomes  significant.  In 
such  cases,  one  of  the  more  rigorous  centroid  MD  algorithms  described  in  the  previous  two 
sections  will  be  necessary. 

In  passing,  we  note  that  three-body  centroid  correlations  (or  ’‘dispersion")  might  be 
incorporated  by  fitting  the  excess  centroid  free  energy  due  to  the  three-body  interaction  to 
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some  functional  form  having  a  few  adjustable  parameters.  This  procedure  is  analogous  to 
the  determination  of  pseudopotentials  from  ab  imtio  quantum  chemistry  calculations. 

V.  NUMERICAL  RESULTS 

In  this  section,  the  accuracy  and  efficiency  of  the  various  centroid  MD  algorithms  is 
examined  for  both  prototype  low-dimensional  models  and  more  realistic  many-body  systems. 
Due  to  limitations  in  the  available  computational  resources,  not  every  model  has  been  studied 
with  every  algorithm.  Nevertheless,  the  numerical  studies  below  clearly  illustrate  the  major 
strengths  and  u’eaknesses  of  each  approach.  AppUcations  of  centroid  MD  to  actual  physical 
systems  will  be  the  subject  of  future  publications. 

The  numerical  examples  and  associated  discussion  are  presented  in  the  order  that  the 
algorithms  appear  in  the  previous  sections.  All  of  the  algorithms  are  essentially  the  same  in 
the  numerical  MD  propagation  of  the  centroid  variable(s),  except  for  the  way  in  which  the 
centroid  force  is  calculated.  As  specified  by  the  centroid  MD  algorithm,  [l,3,4j  the  initial 
centroid  positions  were  generated  by  Metropolis  importance  sampling  from  the  centroid 
distribution,  while  the  initial  momenta  were  directly  sampled  from  the  Gaussian  centroid 
momentum  distribution.  In  most  cases,  the  centroid  correlation  functions  were  plotted 
along  with  the  exact  Kubo  transformed  correlation  function.  [161  if  available.  If  required, 
the  Fourier  relation  in  Eq.  (1.1)  can  be  used  to  convert  the  centroid  correlation  functions  to 
the  usual  ones.  [3] 

A.  Normal  Mode  Path  Integral  Molecular  Dynamics 

In  order  to  test  of  the  convergence  of  the  algorithms  on  a  simple  model  for  which  exact 
(luantum  dynamics  can  be  obtained,  the  same  one-dimensional  anharmonic  potential  has 
been  studied  as  in  Paper  II.  This  potential  is  given  by 

+  9 (S.l) 
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with  the  parameters  c  =  0.10.  g  =  0.01.  h  =  1.0,  and  m  =  1.0.  In  these  units,  the  inverse 
temperature  J  is  given  in  terms  of  the  dimensionless  parameter  3hu.  The  potential  has  a 
single  minimum  at  g  =  0.  exhibiting  a  cubic  anharmonicity  for  small  deviations  from  the 
minimum  and  a  quartic  anharmonicity  for  larger  deviations  from  the  minimum.  At  low 
temperatures  where  mode  quantization  effects  become  important,  the  cubic  anharmonicity 
is  the  dominant  perturbation.  The  anharmonicities  in  Eq.  (5.1)  perturb  the  energy  gap 
between  the  ground  and  first  excited  vibrational  state  downward  by  6%  from  the  harmonic 
limit.  This  corresponds  to  a  large  anharmonicity  for  a  typical  molecular  vibration.  The 
temperature  employed  in  the  calculations  is  d  =  5  which  is  typical  of  a  C— C  single  bond  at 
300  K.  For  this  one-dimensional  model,  the  exact  results  were  obtained  by  first  diagonalizing 
the  nonlinear  Hamiltonian  in  a  harmonic  oscillator  basis  set  to  find  its  eigenvalues  and 
eigenvectors.  Then,  one  hundred  eigenstates  were  employed  to  calculate  the  exact  Kubo- 
transformed  time  correlation  fimctions.  [3] 

In  the  first  case,  the  centroid  force  was  calculated  by  normal  mode  path  integral  MD 
averaging  at  each  time  step.  As  discussed  in  Sec.  II,  the  fictitious  masses  of  all  P  - 1  normal 
modes  except  the  centroid  mode  were  chosen  to  be  rUn  =  rn!/0X'^,  where  m'  =  (3,  so  that 
the  normal  modes  all  oscillate  with  unit  frequency.  Nose- Hoover  chains  consisting  of  two 
thermostats  with  fictitious  masses  of  1.0  and  a  fictitious  temperature  of  1.0  were  attached 
to  each  path  integral  normal  mode.  [24]  The  multiple  time  step  method  [29]  allowed  the 
integration  of  the  motion  due  to  the  residue  potential  AV{q)  for  every  10  steps  of  linear 
motion  of  the  reference  harmonic  path  integral  normal  modes.  The  centroid  motion  was 
integrated  in  terms  of  the  normal  modes  using  the  velocity  Verlet  algorithm  with  a  time 
step  of  0.05.  The  Nose-Hoover  chain  dynamics  were  calculated  for  10“^  steps  with  a  time- 
•step  of  0.005.  and  the  centroid  force  [Eq.  (2.1)],  as  well  as  the  centroid  curvature  [Eq. 
(2.18)]  and  force  fluctuations  [Eq.  (2.24)].  were  accumulated  during  the  simulation.  For  the 
position  correlation  function  of  the  prototype  model  in  Eq.  (5.1),  the  dynamics  of  10'*  centroid 
trajectories  were  averaged  unless  otherwise  specified.  To  speed  up  the  convergence  of  the 
centroid  force,  the  frequency  of  the  harmonic  reference  potential  was  constantly  updated 
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according  to  Eq.  ^2.16).  Moreover,  since  the  force  changes  smoothly  with  the  centroid 
motion,  the  centroid  force  could  be  evaluated  via  the  gradient  correction  in  Eqs.  (2.17)  and 
(2.18)  if  the  change  were  within  some  tolerance.  By  setting  this  tolerance  to  20%-beyond 
which  a  new  .\MPIMD  calculation  of  the  centroid  force  was  initiated— the  speed  of  the 
centroid  MD  integration  could  be  increased  by  a  factor  of  two  to  three.  The  use  of  a  smaller 
fictitious  mass  in  the  spirit  of  an  extended  Lagrangian  simulation  was  not  implemented  in 
the  present  work,  but  such  an  approach  has  been  explored  recently  with  good  success  by 
Martyna  for  centroid  MD.  [26] 

In  Fig.  1.  the  centroid  position  correlation  function  is  shown  for  the  nonlinear  oscillator 
of  Eq.  (5.1)  at  J  =  5.  The  solid  circles  are  the  exact  quantum  results,  the  solid  line  is  the 
NMPIMD  result  for  the  discretization  parameter  of  P  =  20.  the  dashed  line  is  the  NMPIMD 
result  for  P  =  10.  and  the  dot-dashed  line  is  the  result  from  the  effective  harmonic  approxi¬ 
mation  for  the  centroid  force.  [3]  The  improvement  due  to  the  increased  value  of  P  is  evident, 
and  the  brute  force  NMPIMD  calculation  is  superior  to  the  effective  harmonic  approxima¬ 
tion.  However,  the  computation  time  increases  substantially  with  increasing  P.  Obviously, 
calculating  the  centroid  force  exactly  using  NMPIMD  at  each  time  step  provides  the  most 
accurate  dynamical  result,  but  the  method  as  presently  applied  will  become  time  consuming 
for  many-body  systems  having  a  large  number  of  beads.  (The  NMPIMD  calculations  in  Fig. 

1  took  approximately  50  CPU  hours  on  an  IBM  RISC/6000  Model  370  workstation.) 

.4s  an  aside,  we  note  that  the  results  shown  in  Fig.  1  provide  further  evidence  that 
centroid  MD  is  a  method  valid  for  general  nonlinear  potentials  and  is  not  intrinsically  "v-ed- 
ded”  to  the  effective  harmonic  perspective.  (The  latter  perspective  simply  helps  to  justify 
centroid  MD.  T.3])  The  results  in  Fig.  1  show  that  a  centroid  MD  calculation  will  be  more 
accurate  when  the  centroid  force  is  computed  exactly,  as  opposed  to  when  it  is  approximated 
by  the  effective  harmonic  expression.  The  analytical  justification  of  centroid  MD  for  general 
systems  is.  of  course,  the  subject  of  the  companion  paper.  :4] 
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B.  Staging  Path  Integral  Monte  Carlo 


As  discussed  in  Sec.  IIB.  another  strategy  for  computing  the  centroid  force  in  centroid 
MD  is  the  the  staging  path  integral  MC  method.  The  specific  staging  algorithm  introduces 
a  stochastic  d\'namics.  which  in  the  limit  of  N\fc  oc  [cf.  Sec.  IIB)  gives  exactly  the 
centroid  force  and  centroid  MD.  To  test  the  efficiency  and  convergence  of  this  method,  the 
centroid  position  correlation  function  for  the  potential  in  £q.  (5.1)  was  calculated  at  =  o.O 
as  in  Pig.  1.  A  harmonic  reference  propagator  of  frequency  cj  =  1.0  and  P  =  10  was  used 
in  the  staging  path  integral  MC  with  the  number  of  trial  chaiin  quasiparticles  in  the  staging 
algorithm  adjusted  to  yield  an  acceptance  rate  of  50%.  The  staging  centroid  MD  calculations 
converged  to  a  degree  comparable  to  that  shown  in  Fig.  1  for  a  value  of  N\fc  =  10.  taking 
only  a  few  CPU  minutes  on  the  IBM  workstation.  This  calculation  is  therefore  considerably 
faster  than  the  one  performed  with  NMPIMD.  Since  the  staging  PIMC  method  approximates 
the  centroid  MD  asymptotically,  the  averaging  of  the  centroid  force  is  much  less  demanding. 
In  addition,  the  normal  mode  transformation  in  the  NMPIMD  algorithm,  which  consumes  a 
large  portion  of  CPU  time  for  a  large  value  of  P,  is  completely  avoided  in  the  staging  PIMC 
method  for  centroid  MD  as  described  in  Sec.  IIB. 

The  above  result  strongly  suggests  the  feasibility  of  using  centroid  MD  to  simulate  the 
dynamical  time  correlation  functions  of  quantum  particles  in  many-dimensional  systems 
As  a  test  case,  the  self-diffusion  process  of  a  quantum  particle  in  a  classical  solvent  was 
“^studied.  The  solvent  in  this  case  was  a  Lennard-Jones  fluid  with  the  parameters  a  =  2.556 
A.  €  =  10.22  K.  and  the  solute-solvent  pair  potential  had  the  simple  form  [30]  vir)  = 
[B/{C  -i-  r^)  -  l](A/r‘‘),  where  A  =  0.665.  B  =  89099,  C  =  12608  (in  atomic  units).  The 
simulation  consisted  of  512  total  particles  at  a  temperature  T  =  309  K  and  a  reduced 
density  of  p’  =  0.3.  The  solute  was  assumed  to  be  a  quantum  particle  having  the  electron 
mass  TTif,  but  Planck’s  constant  h  was  reduced  by  a  factor  of  10  so  that  a  discretization 
number  of  F  =  100  was  adequate  in  the  simulation.  The  dynamics  of  200  centroid  MD 
trajectories  were  calculated  for  150  time  steps  with  an  increment  of  0.01  with  '^Arc  =  10 
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in  the  staging/ centroid  MD  algorithm.  All  of  the  data  was  scaled  by  the  Lennard- Jones 
potential  width  a.  depth  e.  and  the  electron  mass  m,.  The  mean-squared  displacement  of 
the  quantum  solute  is  plotted  in  Fig.  2  for  the  centroid  MD  with  staging  PIMC  method 
(solid  line)  and  for  classical  MD  (dashed  line).  The  quantum  effect  is  clearly  manifested, 
though  the  quantum  solute  is  far  from  localized  [30,36].  This  calculation  took  approximately 
one  day  of  dedicated  CPU  time  on  the  IBM  RISC/6000  Model  370  workstation. 

Though  quite  promising,  the  staging/centroid  MD  results  are  certainly  preliminary 
Longer  computer  runs  and  better  statistics  are  required  to  evaluate  diffusion  constants  to 
a  acceptable  degree  of  accuracy.  Furthermore,  to  simulate  an  actual  solvated  electron  will 
require  a  \*alue  of  P  =  1000  to  converge  the  quantum  effects.  Nevertheless,  this  exam¬ 
ple  demonstrates  it  is  possible  to  carry  out  centroid  MD  calculations  for  realistic  physical 
systems  using  the  "brute  force’  algorithm  based  on  staging  PIMC. 

C.  Effective  Harmonic  Method:  Extended  Lagrangian  Approach 

The  effective  harmonic  method  for  centroid  MD  [cf.  Sec.  IIIA]  can  be  implemented 
with  the  e.xtended  Lagrangian  approach  of  Sec.  IIIB.  This  approach  will  produce  centroid 
MD  if  the  fictitious  mass  is  sufficiently  small  so  that  the  fictitious  variable,  namely  the 
Gaussian  width  factors  Cc  in  the  present  case,  rapidly  average  the  dynamics  about  the 
minimized  centroid  potential  surface.  As  a  test  of  this  approach,  the  nonlinear  potential  in 
'"Eq.  (5.1)  was  studied  to  determine  the  dependence  of  the  extended  Lagrangian  method  on 
the  fictitious  mass  and  the  time-step  At.  In  Fig.  3,  a  single  centroid  trajectory  is  plotted 
for  TUf  =  0.1  and  At  =  0.01.  Shown  is  the  total  energy  (dot-dashed  line),  the  centroid 
kinetic  energy  (solid  line)  and  the  fictitious  kinetic  energy  on  an  enlarged  scale  by  a  factor 
of  10^  (dotted  line).  Clearly,  under  the  chosen  conditions,  the  fictitious  dynamics  have  little 
effect  on  centroid  dynamics  except  to  generate  the  centroid  mean  force.  For  this  study, 
stable  results  were  obtained  until  the  values  =  0.5  and  At  =  0.1  were  used.  The  "rule 
of  thumb"  is  to  choose  me  small  enough  so  that  the  fictitious  kinetic  energy  is  less  than  one 
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percent  of  the  total  energy  and  to  make  At  small  enough  so  that  the  total  energy  does  not 
drift  more  than  a  few  percent  at  the  end  of  a  single  trajectory.  Although  not  shown  here, 
the  centroid  position  correlation  function  calculated  with  the  extended  Lagrangian  method 
was  found  to  be  the  same  as  the  effective  harmonic  result  shown  in  Fig.  1  (the  latter  result 
was  obtained  by  representing  the  effective  harmonic  centroid  potential  on  a  one-dimensional 
grid).  Since  for  this  example  (and  the  following  one)  the  Gaussian  averages  can  be  expressed 
in  closed  form,  the  extended  Lagrangian  calculation  takes  very  little  CPU  time.  Thus,  the 
numerical  effort  implicit  in  this  method  is  negligible  compared  to  either  of  the  direct  (brute 
force)  path  integral  centroid  MD  methods. 

-As  a  more  challenging  test,  the  effective  harmonic  theory  with  the  extended  Lagrangian 
algorithm  was  applied  to  a  three-dimensional  particle  in  a  non-separable  potential  well,  given 
by 


1  ^ 

HQu  92,^3)  =  +99i)  +  c  919293 

^  i=l 


(5.2) 


with  the  parameters  9  =  0.1,  ft  =  1.0.  m  =  1.0.  and  0  =  5.0.  The  extended  Lagrangian 
centroid  dynamics  were  carried  out  for  the  three  centroid  \’ariables  and  the  six  dynamical 
elements  of  the  Gaussian  width  factor  matrix  with  a  time  step  of  At  =  0.01  and  a  fictitious 
mass  of  me  =  t).l.  The  quantum  position  correlation  function  (not  the  centroid  correlation 
function)  is  shown  in  Fig.  4  for  c  =  0.1  along  with  the  classical  MD  result.  The  e.xtended 
Lagrangian  simulation  for  this  example  was  both  stable  and  efficient,  and  the  quantum 
—  effects  are  seen  to  be  significant. 


D.  Effective  Harmonic  Method:  Iterative  Approach 

The  method  is  rather  straightforward  and  was  the  original  algorithm  used  in  the  previous 
calculations  published  in  Papers  II  and  III.  In  fact,  the  iteration  method  works  extremely 
well  for  the  examples  discussed  in  the  previous  subsection  ( i.e..  even  better  than  the  extended 
Lagrangian)  and  uill  not  be  reproduced  here.  Instead,  the  focus  will  shift  to  a  more  realistic 
many-body  system. 
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As  has  been  pointed  out  in  Sec.  IIIB.  the  most  serious  potential  difficulty  associated 
with  solving  the  transcendental  equations  in  Eqs.  (3.1)— (3.3)  in  many-dimensional  space  is 
the  evaluation  oi  the  Gaussian  averages.  Several  approaches  were  proposed  and  formulated 
in  Sec.  III.  For  the  present  purposes,  the  Gaussian  representation  of  pair  potentials  in  Eq. 
(3.17)  is  the  most  convenient  and  efficient.  A  set  of  five  of  Gaussian  parameters  was  found 
to  be  accurate  in  representing  a  Lennard-Jones  pair  potential.  [37]  In  order  to  demonstrate 
the  validity  of  the  effective  harmonic  approximation  and  the  Gaussian  representation  of  the 
Lennard-Jones  potential,  in  Fig.  5  the  centroid  force  thus  calculated  is  plotted  (solid  line) 
along  with  the  results  of  an  exact  PIMC  calculation  (solid  circles).  The  parameters  were 
taken  to  be  d  =  1.0.  m  =  1.0.  and  h  =  0.1.  so  the  PIMC  calculations  were  performed  for 
P  =  20  with  10^  MC  passes.  The  results  are  in  remarkably  good  agreement. 

Next  the  effective  harmonic  theory  with  the  iterative  algorithm  for  the  centroid  force 
was  applied  to  study  a  quantum  particle  solvated  in  a  pairwise  Lennard-Jones  classical  fluid 
similar  to  that  described  in  Sec.  VB.  The  Lennard-Jones  interaction  between  the  quantum 
particle  and  the  solvent  particles  was  also  reprinted  by  five  independent  Gaussian  func¬ 
tions.  The  well  depth  and  width  of  Lennard-Jones  potential  were  taken  as  the  natural  units, 
and  a  unit  mass  was  assumed  for  both  the  solvent  and  the  solute.  X  total  of  125  solvent 
particles  in  a  periodic  cubic  were  used  to  model  the  solvent  at  a  reduced  density  0.8  and 
at  a  reduced  temperature  1.0  The  quantum  nature  of  the  solute  particle  was  adjusted  by 
varying  the  Planck's  constant  h  from  A  =  0.0  to  ^  =  0.3  in  the  reduced  units.  The  dynamics 
of  10^  independent  centroid  trajectories  were  then  integrated  for  500  time  steps  of  duration 
0.01.  At  the  initial  stage  of  the  centroid  MD  simulation,  10  —  20  iterations  were  needed  for 
convergence  in  computing  the  centroid  force  at  each  time  step.  After  the  centroid  began  to 
move,  the  self-consistent  parameter  set  from  the  previous  time  step  was  taken  as  the  input 
to  the  iterations  at  the  current  time  step  and.  accordingly,  only  1  -  3  iterations  were  needed 
to  converge  the  centroid  force. 

The  mean  square  displacement  of  the  solute  is  plotted  in  Fig.  6  for  centroid  MD  (solid 
line)  and  for  classical  MD  (dashed  line).  The  diffusion  constant  measured  from  the  slope 
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of  the  displacement  c\irve  shows  a  20%  decrease  due  to  quantum  effect  in  this  case.  The 
effecti\'e  harmonic  centroid  MD  method  with  the  iterative  approach  was  found  to  be  stable 
iind  highly  efficient.  The  results  in  Fig.  b  were  obtained  in  several  hours  of  CPU  time  on 
the  IBM  RISC/' 6000  Model  370  workstation. 

E.  Centroid  Pseudopotential  Method 

First  of  all.  it  is  important  to  examine  the  suggestion  in  Sec.  IV  that  an  effective  pairwise 
centroid  potential  can  capture  the  dominant  quantum  effects  in  the  many-body  centroid 
potential  surface.  Since  only  the  force  is  relevant  to  centroid  MD.  the  following  results  will 
be  expressed  in  terms  of  the  centroid  force  instead  of  the  centroid  potential.  For  convenience, 
a  Lennard- Jones  potential  was  employed  with  the  reduced  parameters  a  =  1.0,  e  =  1.0, 
m  =  1.0,  0  =  1.0,  and  h  =  1.0.  The  computational  procedure  is  described  as  follows: 

(1)  The  centroid  force  between  a  pair  of  quantum  particles  interacting  via  the  Lennard 
Jones  potential  was  calculated  as  a  function  of  separation  between  their  centroids  in  a 
PIMC  simulation. 

(2)  The  total  centroid  force  between  the  three  particles  arranged  on  ein  equilateral  triangle 
was  then  calculated  in  a  similar  fashion  by  uniformly  changing  the  distance  of  the  three 
sides  of  the  triangle. 

(3)  The  pairwise  pseudopotential  method  was  then  applied  to  the  same  three-particle  system 
-  with  the  pair  potential  obtained  from  step  (1). 

The  results  from  the  direct  simulation  in  step  (2)  and  the  pseudopotential  approximation 
in  step  (3)  are  plotted  in  Fig.  7  and  agree  quite  well  except  at  very  close  (and  improbable) 
distances. 

In  paper  II.  the  self-diffusion  process  in  fluid  neon  [38]  was  studied  with  centroid  .MD 
by  including  the  leading  correction  in  the  pair  centroid  potential.  The  same  calculation  has 
been  carried  out  here  using  the  pseudopotential  method.  The  parameters  of  the  system  and 
details  of  the  centroid  MD  simulation  can  be  found  in  Paper  II  let  Sec.  IV  C  and  Table  I 
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in  that  paper).  To  illustrate  the  difference  between  the  two  approximations,  in  Fig.  8  the 
pair  centroid  pseudopotential  force  used  in  this  calculation  is  plotted  (solid  line)  against  the 
approximate  centroid  force  used  in  the  Paper  II  (dot-dashed  line)  and  the  classical  Lennard- 
Jones  force  (dashed  line).  Then,  in  Fig.  9  the  centroid  velocity  time  correlation  function 
is  plotted  for  quantum  neon  (solid  line)  using  the  pseudopotential  method  and  for  classical 
neon  (dashed  line).  Evidently,  the  more  repulsive  centroid  force  shown  in  Fig.  8  leads  to 
the  decrease  in  the  quantum  diffusion  constant  seen  in  Fig.  9  (i.e..  through  the  integral  of 
the  centroid  velocity  autocorrelation  function  in  the  Green-Kubo  formula  for  the  diffusion 
constant  [3]). 


VI.  CONCLUDING  REMARKS 

This  paper  is  solely  devoted  to  the  development  of  numerical  algorithms  for  centroid  MD 
and.  in  turn,  it  sets  the  stage  for  a  broad  spectrum  of  applications.  Though  centroid  MD 
makes  the  study  of  quantum  dynamics  in  raany-body  systems  a  possibility,  the  numerical 
integration  of  the  centroid  MD  equations  is  by  no  means  trivial.  This  is  because  the  centroid 
force  Fciqc)  in  Eq,  (1.4)  is  a  quantum  mean  force  which  caimot  be  evaluated  exactly  by  simple 
analytical  expressions.  Therefore,  several  algorithms  for  centroid  MD  have  been  outlined  and 
tested  in  the  present  paper.  Though  these  algorithms  have  been  described  in  considerable 
detail,  some  remarks  are  in  order  regarding  the  outlook  for  their  application  in  realistic 
simulations. 

Strictly  speaking,  the  quantum  imaginary  time  path  fluctuations  about  the  centroid  vari¬ 
able  must  be  completely  averaged  in  order  to  compute  the  centroid  force  in  Eq.(1.4)  at  each 
time  step  in  the  centroid  MD  integration.  The  NMPIMD/centroid  MD  algorithm  in  Sec. 
IIA  will  essentially  produce  the  e.xact  centroid  force  at  each  time  step,  but  the  prelimi¬ 
nary  results  reported  in  Sec.  VA  suggest  its  application  to  a  real  many-body  system  will  be 
formidable  in  any  practical  sense.  Perhaps  a  more  complete  e.xamination  of  the  extended 
Lagrangian  aspects  of  the  NMPIMD  algorithm  will  provide  a  more  optimistic  prognosis  for 
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the  method.  [261  On  the  other  hand,  the  results  from  the  staging  PIMC/centroid  .MD  al¬ 
gorithm  developed  in  Sec.  IIB  suggest  that  a  brute  force  centroid  MD  calculation  is  indeed 
within  the  realm  of  possibility  for  realistic  systems.  The  staging  PIMC/centroid  MD  scheme 
is  simpler  and  faster  than  NMPIMD /centroid  MD  algorithm.  The  preliminary  results  sug¬ 
gest  that  a  single  quantum  particle  immersed  in  a  classiced  solvent  with  P  =  100  -  1000.  or 
an  ensemble  of  quantum  particles,  each  with  F  20  —  30.  can  be  treated  with  the  staging 
PIMC/centroid  MD  algorithm  on  a  higher-end  computer  workstation.  Of  course,  computer 
performance  is  rapidly  increasing,  so  optimism  is  clearly  in  order.  Moreover,  the  outlook  for 
implementing  centroid  MD  on  a  massively  parallel  computer  seems  most  promising,  though 
such  algorithms  remain  to  be  explored. 

For  systems  which  are  presently  beyond  the  reach  of  the  brute  force  path  integral  cen¬ 
troid  MD  methods,  the  effective  harmonic  centroid  MD  algorithm  of  Sec.  Ill  is  an  attractive 
option.  In  fact,  this  algorithm  will  always  provide  a  highly  efficient  alternative  to  the  brute 
force  methods  in  systems  where  the  quantum  modes  are  predominantly  influenced  by  pos¬ 
itive  curvature  in  the  potential  (i.e..  so  that  the  minimization  of  the  effective  harmonic 
potential  is  efficient).  The  vibrational  dynamics  and  relaxation  of  condensed  phase  poly¬ 
atomic  molecules,  or  perhaps  the  dynamics  of  quantum  clusters,  may  be  ideal  candidates 
for  the  effective  harmonic  algorithm. 

Due  to  its  simplicity  and  flexibility,  the  centroid  pseudopotential  method  described  in 
Sec.  IV  is  panicuiarly  attractive  for  very  large  systems,  provided  those  systems  are  not  too 
quantum  in  nature  or  have  low  densities.  For  example,  the  transport  coefficients  of  nearly 
classical  fluids,  simulations  of  quantum  friction,  and  calculations  of  quantum  solvent  spectra 
could  all  be  investigated  with  this  method  (and  probably  with  the  others  as  well— at  least 
in  time). 

As  a  final  point,  we  note  that  the  algorithms  described  in  the  present  paper  might  be 
combined  to  greatly  enhance  the  capability  of  centroid  MD  for  the  study  of  real  systems. 
Consider  the  example  of  proton  exchange  dynamics  between  two  solutes  in  a  quantized 
solvent  such  as  water.  The  proton  centroid  dynamics  might  be  treated  by  one  of  the  brute 
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force  path  integral  methods,  the  solute  intramolecular  centroids  could  be  integrated  through 
the  effective  harmonic  approach,  while  the  solvent-solute  and  solvent-solvent  interactions 
could  be  approximated  by  the  centroid  pseudopotential  method.  The  teasibility  of  such  an 
algorithm  is  insured  by  the  so-called  free  energy  Born-Oppenheimer  approximation.  (19)  A 
\'ariety  of  applications  of  centroid  MD  will  be  reported  in  future  publications. 

ACKNOWLEDGMENTS 

This  research  was  supported  by  grants  from  the  National  Science  Foundation  and  the 
Office  of  Naval  Research.  GAV  is  a  recipient  of  a  National  Science  Foundation  Presidential 
Young  Investigator  Award,  a  David  and  Lucile  Packard  Fellowship  in  Science  and  Engineer¬ 
ing,  an  Alfred  P.  Sloan  Foundation  Research  Fellowship,  and  a  Dreyfus  Foundation  New 
Faculty  Award.  The  authors  are  indebted  to  Diane  Sagnella  for  her  critical  reading  of  this 
manuscript. 


37 


APPENDIX  A;  EVALUATION  OF  THE  SUBMATRIX  IN  EQ.  (3.13) 


In  general,  the  centroid-constrained  correlation  function  in  Eq.  (3.3)  and  Eq.  (3.5)  is  a 
{3N  X  3A’')— th  order  matrix  for  an  A’— body  system  in  three  dimensional  space.  For  the 
simplicity  of  notation,  the  three  space  indices  will  be  eliminated  here.  Thereby,  for  a  particle 
with  three  spatial  components  the  matrix  Cc  is  re-defined  below  as  an  iV— th  order  matrix, 
with  each  element  understood  to  be  defined  as 

^  ni2  /--i3  \ 

^cAj 

c„,  3  C?‘,  CSi  ■  (Al) 

^31  /^32  /^33 

y  ^c.ij  j 

Here,  the  subscripts  denote  the  particle  index  and  the  superscripts  denote  the  three  spatial 
dimensions.  It  is  also  to  be  noted  that  all  elements  of  this  matrix  depend  on  the  centroid 
variables  of  the  iV— body  system.  One  can  now  partition  the  inverse  of  the  matrix  Cc, 
defined  as  B,  such  that 


=  B  = 


bi2 

^21  t>22 


Here,  bn  =  Bu  is  treated  as  a  matrix  element  [which  is.  according  to  Eq.  (Al),  actually  a 
matrix  of  order  3l.  =  (B12,  B13,  ■  ■■ .  Buv)  is  like  a  row  vector.  b2i  is  like  a  column  vector 

such  that  boi  =  bi2.  and  b22  is  like  an  (iV  -  1)  x  (iV  -  l)-th  order  matrix. 

With  the  above  notation  in  hand,  all  of  the  other  variables  {qrj,  —  qj^}  in  Eq.  (3.8) 
can  be  integrated  out  to  give  a  reduced  Gaussian  average  [cf.  Eq.  (3.12)],  with  the  reduced 
Gaussian  distribution  given  by 

P(o  )  =  exp(-q-B-q/2) 

*  fdqexp(-q-B-q/2) 

det(27rbooh  1  -  ,  ,  .  _i,  . 

=  expi- 7i(6u  -  bi2bo2b2U?i/2]  AS) 

On  the  other  hand,  one  has  the  matrix  element  from  Eq.  (.A.1),  i.e.. 


Cr.ii  =  (B“‘)ii  =  det(b22j/ det(B) 
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where  the  denominator,  after  some  algebra,  can  be  written  as 


det(B)  =  (611  —  bi2b22^b2i)det(bo2) 


CAS) 


Therefore,  the  reduced  matrix  C'  =  Ccai  can  be  introduced  so  that  the  reduced  Gaussian 
distribution  in  Eq.  (A3)  reads  as 


“  diJ(dcp)  9-  ■  •  «i/2)  ■  (A6) 

where  according  to  the  notation  used  here  C'  is  a  matrix  of  order  3  and  qi  is  a  three- 
dimensional  vector. 

.\s  stated  in  the  main  text,  a  pair  potential  only  depends  on  the  relative  position  of  the 
pair  particles.  Therefore,  the  preceding  result  should  be  expressed  in  terms  of  =  qi—qj, 
the  vector  connecting  particles  i  and  j.  Thereby,  an  uniform  transformation  is  required  to 
obtain  the  relevant  matrix  element,  i.e.. 


p-l  _  9f>-l  I 


(A7) 


which  defines  the  reduced  Gaussian  distribution  for  the  variable  in  Eqs.  (3.11)  and  (3.12). 
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FIGURES 

FIG.  1.  A  plot  of  the  centroid  position  correlation  function  for  the  nonlinear  oscillator  of  Eq. 
(5.1)  at  ii  =  5.  The  solid  circles  are  the  exact  quantum  results,  the  solid  line  is  the  NMPIMD 
result  for  P  =  20.  the  dashed  line  is  the  NMPIMD  result  for  P  =  10.  and  the  dot-dashed  line  is 
the  result  of  the  effective  harmonic  approximation  for  the  centroid  force  (see  Paper  II) 

FIG.  2.  A  plot  of  the  mean  square  displacement  for  a  quanttun  particle  solvated  in  a  classical 
neon  fluid.  The  solid  line  is  the  stochastic  centroid  MD  with  staging  PIMC  described  in  Sec.  IIB 
while  the  dashed  line  is  for  classical  MD. 

FIG.  3.  A  plot  of  the  energy  data  for  a  single  trajectory  for  the  potential  in  Eq.  (5.1)  as 
calculated  with  the  extended  Lagrangian  centroid  MD  method  of  Sec.  IIIB.  The  dot-dashed  line  is 
the  total  energy,  the  solid  line  is  the  centroid  kinetic  energy,  and  the  dashed  line  is  the  fictitious 
kinetic  energy  enlarged  by  a  factor  of  10^. 

FIG.  4.  A  plot  of  the  real  time  position  correlation  function  calculated  by  the  extended  La¬ 
grangian  centroid  MO  method  of  Sec.  IIIB  for  the  three-dimensional  potential  given  in  Eq.  (5.2) 
at  a  temperattire  oi  3  =  5.  Also  shown  is  the  classical  MD  result. 

FIG.  5.  plot  of  the  centroid  force  for  a  Lennard-Jones  potential.  The  solid  circles  represent 
e.xact  results  from  a  PIMC  simulation,  while  the  solid  line  represents  the  effective  harmonic  approx¬ 
imation  from  Eqs.  ^^3.1)— (3.3)  with  the  Gaussian  representation  of  the  Lennard-Jones  potential  in 
"Eq.  (3.17). 

FIG.  6.  .A.  plot  of  the  mean  squared  displacement  of  a  quantum  Lennsu'd-Jones  particle  solvated 
in  a  classical  Lennard-Jones  fluid.  The  solid  line  is  the  centroid  MD  result  obtained  with  the 
effective  harmonic  centroid  force  solved  iteratively  (cf.  Sec.  IIIC).  while  the  dashed  line  is  the 
classical  MD  result. 
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FIG.  7.  A  pi  01  of  the  total  centroid  force  between  three  Lennard- Jones  particles  distributed  on 
a  equilaterad  triangle.  The  solid  circles  are  the  exact  PIMC  simulation  results,  while  the  solid  line 
is  the  pseudopotential  approximation  of  Sec.  IV. 

FIG.  8.  A  plot  of  the  pair  force  between  two  Lennard- Jones  particles  as  a  function  of  separation. 
The  solid  line  is  exact  pair  centroid  force  used  in  the  pseudopotential  calculation,  the  dot-dashed 
line  is  the  approximate  pair  centroid  force  used  in  Paper  II.  and  the  dashed  line  is  the  classical 
Lennard-Jones  force. 

FIG.  9.  A  plot  of  the  centroid  velocity  autocorrelation  function  for  fluid  neon.  The  solid  line 
is  the  centroid  MO  result  calculated  with  the  centroid  pseudopotential  approximation,  while  the 
dashed  line  is  the  classical  MD  result. 
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of  a  Solvated  Quantum  Particle 
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